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ABSTRACT 

The  modeling  of  damped  signals  as  the  impulse  response  of  a  pole- zero  system 
is  considered  for  a  broad  range  of  pole-zero  modeling  algorithms.  The  goal  is  to 
obtain  the  best  possible  fit  between  the  model  impulse  response  and  the  modeled  sig- 
nal. Prony's  method,  the  least  squares  modified  Yule- Walker  equations  (LSMYWE), 
iterative  prefiltering,  and  the  Akakie  maximum  likelihood  estimator  are  compared 
on  known  test  sequences  for  a  variety  of  model  degrading  situations  (e.g.,  additive 
noise)  to  develop  an  understanding  of  which  methods  are  most  suitable  for  mod- 
eling real  world  signals.  A  correlation  domain  version  of  interative  prefiltering  is 
also  introduced.  The  most  robust  algorithms  are  determined  to  be  LSMYWE  using 
singular  value  decomposition  and  iterative  prefiltering  (including  the  correlation  do- 
main version).  Modeling  several  laboratory  generated  short  duration  acoustic  signals 
confirmed  the  robustness  of  LSMYWE  and  iterative  prefiltering.  It  is  shown  that 
correlation  domain  iterative  prefiltering  outperforms  standard  iterative  prefiltering 
when  large  model  orders  are  required  for  accurate  modeling.  Shank's  method  was 
determined  to  be  the  most  effective  method  of  determining  the  zeros  of  a  pole-zero 
model  when  a  time  domain  match  is  required. 
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I.  INTRODUCTION 

A.      RATIONAL  MODELING  OF  TIME  SERIES  DATA 

The  need  to  find  a  compact  parametric  representation  for  time  series  data  arises 
in  many  fields  of  study.  In  electrical  engineering,  finding  such  representations,  or 
models,  appears  under  such  topics  as  optimum  control,  optimum  filtering,  system 
identification,  model  order  reduction,  waveform  encoding,  and  spectrum  estimation. 
Other  fields  refer  to  such  modeling  as  time  series  analysis  or  forecasting.  In  digital 
systems,  the  model  parameters  take  the  form  of  difference  equation  coefficients  or, 
alternatively,  transfer  function  coefficients.  The  most  general  form  of  a  difference 
equation  is  one  that  uses  both  feed-forward  and  feedback  coefficients.  The  corre- 
sponding transfer  function  is  a  ratio  of  two  polynomials  in  the  complex  variable  z 
and  is  refered  to  as  a  rational,  or  pole-zero,  model. 

Historically,  the  more  general  pole-zero  model  has  been  used  in  relatively  few 
applications  compared  to  all-pole  (feedback  only)  and  nil-zero  (feed  forward  only) 
models.  In  some  cases,  an  all-pole  or  all-zero  model  is  the  most  appropriate.  More 
often,  however,  the  all-pole  or  all-zero  model  is  chosen  because  the  optimal  model 
estimation  procedures  are  better  understood  and  easier  to  implement  than  those  for 
pole-zero  modeling.  This  is  particularly  true  for  the  all-pole  case  which,  in  many 
situations,  can  be  determined  by  solving  a  set  of  linear  equations.  Two  recent  devel- 
opments have  led  to  increased  activity  in  applying  pole-zero  models:  (1)  technological 
advances  in  digital  hardware  have  dramatically  reduced  computation  costs  and,  (2) 
a  greater  variety  of  efficient  techniques  for  estimating  pole-zero  model  parameters  is 
now  available. 


Literally  hundreds  of  pole-zero  modeling  estimation  algorithms  and  applications 
have  been  published.  The  majority  of  these  are  built  around  probablistic  or  stochas- 
tic modeling  techniques.  This  is  because  stochastic  modeling  is  usually  the  most 
appropriate  to  forecasting  [Ref.  1,  2,  3]  and  spectrum  estimation  [Ref.  4,  5,  6]  where 
very  little  is  known  about  the  system  input  which  produced  the  time  series  being 
modeled.  A  deterministic  methodology  has  also  been  used  in  which  the  input  and 
output  time  series  are  available  and  the  linear  system  which  'best'  (usually  in  a  least 
squares  sense)  produces  this  cause  and  effect  is  determined.  This  approach  is  usually 
found  under  the  topics  system  identification  [Ref.  7,  8,  9,  10]  and  waveform  encoding 
[Ref.  11].  Another  large  body  of  literature  which  exists  in  parallel  to  stochastic  and 
deterministic  modeling  is  that  of  reduced-order  modeling  [Ref.  12]  which  largely  deals 
with  the  system  control  applications  of  pole-zero  modeling. 

Although  stochastic  and  deterministic  pole-zero  modeling  have  the  same  goal 
and  use  the  same  mathematical  techniques,  the  distinction  is  significant  in  the  way 
data  is  treated  and  in  the  performance  criterion  postulated.  Specifically,  stationarity 
requirements  and  assumptions  about  the  probability  density  functions  of  random  pro- 
cesses in  stochastic  modeling  can  be  limiting  when  dealing  with  real  world  systems.  In 
contrast,  deterministic  techniques  make  no  specific  assumptions  about  the  time  series 
to  be  modeled  except,  of  course,  that  the  form  of  the  model  chosen  is  appropriate. 

One  particular  deterministic  modeling  problem  that  has  received  little  detailed 
attention  is  that  of  finding  a  pole-zero  model  when  the  time  series  being  modeled 
is  a  transient,  'impulse  response'-like  waveform.  Examples  of  situations  where  such 
models  may  be  useful  include  in  modal  or  shock  analysis  of  mechanical  systems  [Ref. 
13,  14],  antenna  response  to  electromagnetic  pulses  [Ref.  15],  and  wavelet  estima- 
tion in  siesmic  studies  [Ref.  16].  Pole-zero  models  make  sense  for  transient  wave- 
forms because  such  a  model  is  structural  for  many  transient  signals.  That  is,  impulse 


response  like  transients  are  usually  the  result  of  a  system  of  oscillators  which  have 
been  excited  for  a  very  short  time  period  relative  to  the  natural  frequencies  of  the 
oscillators.  Pole  pairs  in  pole-zero  models  correspondingly  represent  the  resonant  fre- 
quencies of  a  digital  system.  System  zeros  allow  the  initial  conditions  or  phasing  of 
the  resonant  frequencies  to  be  modeled.  This  combination  of  poles  and  zeros  allows 
signals  to  be  modeled  based  on  time  domain  matching.  When  an  effective  time  do- 
main match  is  achieved  many  concerns  about  assumptions  during  modeling  become 
moot;  an  effective  time  domain  match  has,  by  definition,  effectively  characterized  the 
signal  in  question.  Previous  work  has  concentrated  on  finding  only  transient  model 
poles  [Ref.  17,  18,  19,  20]  or  concentrated  on  one  particular  technique  of  pole-zero 
impulse  response  matching  [Ref.  21]. 

B.     THESIS  OUTLINE 

This  thesis  provides  a  performance  comparison  of  several  pole-zero  modeling 
procedures  applied  to  the  problem  of  model  estimation  from  impulse  response  data. 
The  procedures  compared  are  chosen  to  form  a  cross  section  of  optimality  and  compu- 
tational complexity  of  available  techniques.  In  Chapter  II,  the  modeling  procedures 
selected  for  study  are  described  in  detail.  The  remaining  three  chapters  are  concerned 
with  comparing  the  performance  of  these  modeling  techniques  and  are  organized  as 
follows: 

1.  To  study  the  specific  modeling  properties  of  each  method,  test  impulse  response 
sequences  are  modeled  in  Chapter  Three.  Test  sequences  are  constructed  to 
simulate  the  degradations  likely  to  be  present  in  'real  world'  transient  signals 
(e.g.  noise,  unknown  model  order). 

2.  Using  the  results  obtained  in  Chapter  III,  laboratory  generated  acoustic  tran- 
sient data  is  modeled  in  Chapter  IV.  This  data  is  considered  to  determine  the 


performance  of  various  techniques  when  modeling  the  highly  complex  data  char- 
acteristic of  real  world  sources  about  which  very  little  is  usually  known. 

3.  Chapter  V  summarizes  the  main  conclusions  drawn  from  the  results  in  Chapters 
III  and  IV.  Recommendations  for  further  study  are  also  presented. 


II.  POLE-ZERO  MODELING 

A.      OVERVIEW— THE  VARIETY  OF  TECHNIQUES 

Choosing  a  pole-zero  modeling  technique  from  among  the  many  available  tech- 
niques can  be  difficult.  The  complexity,  applicability,  and  demonstrated  effectiveness 
of  the  different  methods  are  not  always  well  documented.  Additional  consideration  of 
the  many  refinements  that  often  evolve  as  a  technique  is  applied  to  different  problems 
can  lead  to  a  perplexing  array  of  tools  with  which  to  attack  the  modeling  problem. 
For  the  particular  case  of  modeling  by  impulse  response  matching,  little  work  has  ap- 
peared which  sorts  out  the  strengths  and  weaknesses  of  available  modeling  techniques 
or,  in  fact,  indicates  which  methods  may  or  may  not  be  applicable. 

The  basic  scheme  for  fitting  a  pole- zero  system  impulse  response  to  a  given  data 
sequence,  x(n),  is  illustrated  in  Figure  2.1a.  This  is  sometimes  referred  to  as  the  direct 
model.  As  we  shall  see,  the  formulation  of  this  problem  leads  to  a  set  of  nonlinear 
equations  which  require  the  use  of  iterative  techniques  to  solve.  To  overcome  the 
complexities  inherent  in  solving  nonlinear  equations,  the  impulse  response  matching 
problem  can  be  reformulated  as  shown  in  Figure  2.1b.  This  may  be  referred  to  as  the 
indirect  method.  Note  that  these  two  formulations  are  not  equivalent:  the  error  of 
the  direct  method  of  Figure  2.1a  is  Cd(n)  =  x(n)  —  h(n)  while  the  error  for  the  indirect 
method  of  Figure  2.1b  is  e,(n)  =  b(n)  —  a(n)*x(n),  where  h(n)  is  the  impulse  response 
of  the  pole-zero  system  being  found,  b(n)  is  the  corresponding  sequence  of  numerator 
coefficients  and  a(n)  is  the  corresponding  sequence  of  denominator  coefficients.  The 
solution  of  the  indirect  problem  is  considered  suboptimal  in  the  sense  that,  except 
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Figure  2.1:  Two  pole-zero  modeling  problem  formulations:  (a)  direct  for- 
mulation and  (b)  indirect  formulation. 


when  the  modeling  error  goes  to  zero,  the  effect  of  minimizing  e,(n)  is  not  the  same 
as  minimizing  ed{n)  in  the  direct  method. 

One  way  to  organize  pole-zero  modeling  techniques  is  shown  in  Figure  2.2.  In 
deterministic  waveform  matching,  the  ideal  equations  relating  the  proposed  model  to 
the  available  input  and  output  data  sequences  are  formed.  The  solution  which  'best' 
satisfies  these  equations  is  chosen.  In  this  context,  'best'  usually  means  the  solution 
which  minimizes  the  sum  of  squares  of  the  equation  error.  These  are  essentially  the 
Prony  type  methods  [Ref.  17,  22]  (when  the  input  is  an  impulse  response)  and  least 
squares  system  identification  [Ref.  7,  8,  23]  methods  (for  general  input  sequences). 
A  number  of  iterative  techniques  for  waveform  matching  have  also  been  proposed. 
Waveform  fitting  error  [Ref.  24,  25]  and  inverse  filtering  error  [Ref.  26,  27]  are  the 
criteria  most  used. 

Linear  stochastic  pole-zero  modeling  techniques  rely  primarily  on  estimates  of 
second  order  statistics  (auto-  and  cross-correlations)  to  estimate  model  parameters. 
Spectral  estimation  has  been  a  driving  force  for  these  methods  which  are  based  on 
solving  some  form  of  the  modified  Yule- Walker  equations  (see  Chapter  III).  Other 
methods  which  utilize  reflection  coefficients  [Ref.  28]  and  higher  order  statistics  [Ref. 
29,  30]  have  also  appeared. 

The  maximum  likelihood  technique  seeks  parameter  estimates  for  which  the  ob- 
served data  is  the  most  probable  in  the  sense  that  its  conditional  probability  density 
function  (likelihood  function)  is  maximized.  This  technique  is  considered  to  be  sta- 
tistically optimum  but  is  quite  difficult  to  use  because  its  implementation  generally 
requires  the  minimization  of  a  highly  nonlinear  function  [Ref.  1,  31,  32,  33].  Four 
widely  applied  methods  from  each  of  the  categories  of  Figure  2.2  will  be  used  in 
this  thesis.  A  number  of  improvements  which  have  subsequently  been  suggested  for 
these  methods  will  also  be  considered.  Most  recent  work  in  pole-zero  modeling  and 


spectrum  estimation  has  occured  at  the  internal  boundaries  of  Figure  2.2,  i.e.  equiv- 
alent linear  techniques  are  sought  that  perform  as  well  as  modeling  formulations  that 
require  solving  nonlinear  equations  [Ref.  34,  35,  28,  36,  37].  The  application  of 
these  newer  techniques  to  transient  modeling  will  not  be  considered  here  since  we 
expect  that  the  performance  of  the  methods  chosen  will  in  most  cases  bracket  the 
performance  of  these  newer  techniques  with  the  possible  sacrifice  of  computational 
efficiency.  The  four  pole-zero  modeling  techniques  chosen  will  be  described  in  the  re- 
mainder of  this  chapter.  The  transient  modeling  performance  of  these  methods  that 
will  then  be  compared  in  Chapters  III  and  IV. 


Linear 


Iterative 


Deterministic 


Equation  Error 
Methods 

Waveform  Matching 
Inverse  Filtering 

Correlation  Equation 
Error  Methods 

Maximum 
Likelihood  Methods 

Stochastic 


Figure  2.2:    One  way  of  organizing  the  various  pole-zero  modeling  tech- 
niques. 


B.     THESIS  MODELING  TECHNIQUES 
1.     Prony's  Method 

One  of  the  best  known  indirect  techniques  for  matching  a  waveform  to 
the  impulse  response  of  a  linear  time-invariant  system  is  Prony's  method  [Ref.  22]. 

8 


This  method  is  in  fact  a  special  case  of  least  squares  system  identification  in  which 
the  system  input  sequence  is  taken  to  be  a  unit  impulse  and  the  numerator  and 
denominator  coefficients  are  determined  separately. 

The  pole-zero  modeling  problem  is  formulated  as  follows.  The  time  domain 
difference  equation  for  a  general  feedback,  feed-forward  system  can  be  written 


P  Q 

^2  ajx{n  -  j)  =  5^  biu(n  -  i) 

j=0  t=0 


(2.1) 


where  u(n)  is  the  input  sequence  and  x(n)  is  the  output  sequence.  When  the  input 
sequence  is  taken  to  be  a  unit  impulse  (unit  sample  function)  and  the  output  is  taken 
to  be  the  corresponding  impulse  response,  (2.1)  can  be  expressed  in  the  form  of  a 
matrix  equation, 


x(0) 
x(l) 


0 
*(0) 


x(N-l)    x(N-2) 


'  b0  ' 

0 

"    1    " 

0 

a\ 

= 

bQ 
0 

x{N-P)  . 

.  aP  . 

0 

(2.2) 


where  N  is  the  number  of  data  points  used  and,  without  loss  of  generalization,  ao  is 
set  equal  to  one. 

Equation  (2.2)  can  be  solved  by  partitioning, 


a  = 


b 
0 


(2.3) 


where  a  =  [1  a^  •  •  •  ap]T,  b  =  [b0  b\  •  •  •  6q]t  and  X^  and  X#  are  the  corresponding 
lower  and  upper  partitions  of  the  data  matrix  in  (2.2).  The  upper  partition  consists 
of  the  first  Q  +  1  rows  of  the  data  matrix  and  the  lower  partition  is  composed  of  the 
remaining  rows. 

The  solution  can  then  be  obtained  by  first  solving  the  lower  partition, 


X,ia  =  0, 
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(2.4) 


for  a  and  then  finding  b  from  the  upper  partition, 

XBa  =  b.  (2.5) 

lfN  =  P  +  Q  +  l  then  (2.2)  may  have  a  unique  solution  and  the  model  impulse 
response  will  exactly  match  x(n)  forO<n<P  +  Q  +  l.  This  is  referred  to  as  the 
Pade  approximation  [Ref.  22]. 

In  most  circumstances,  however,  the  length  of  the  available  data  sequence 
far  exceeds  P  +  Q  +  1.  It  is  then  desirable  to  use  all  available  data  in  setting  up  (2.2). 
This  leads  to  an  overdetermined  set  of  linear  equations  for  the  lower  partition.  No 
exact  solution  to  (2.4)  usually  exists  in  this  case.  The  relationship  in  (2.4)  becomes 

X^a  =  e,  (2.6) 

where  e  is  the  equation  error  that  will  be  present.  The  solution  of  (2.4)  and  (2.6) 
requires  the  partitioning  of  X,4  as  follows, 

*a=[xa    \    X'A]  (2-7) 

where  x^  is  the  first  column  of  X^.  If  the  remaining  matrix  X^  is  square  and  of  full 
rank,  then  the  solution  to  (2.4)  is  given  by 

a'  =  (XUr'xA  (2.8) 

where 


a=     a,     .  (2.9) 

Otherwise  the  least  squares  error  solution  of  (2.6),  which  minimizes  the  squared  error 
eTe,  is  given  by 

a'  =  X'+xA  (2.10) 
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where  X.'^   is  the  Moore- Penrose  pseudoinverse  of  X.'A.  When  X^  is  of  full  (column) 
rank  P,  the  psuedoinverse  is  given  by 


a'  =  (X'JX'J^XW  (2.11) 

(See  [Ref.  4,  pp.  28-33]  for  an  example  of  this  derivation).  Otherwise,  the  pseudoin- 
verse is  defined  by 

w     t 

where  <7t,  u,,  and  v,  are  defined  by  the  singular  value  decomposition  (SVD)  represen- 
tation of  X^, 

w 

X^  =  X>,u,v,T.  (2.13) 

t=0 

The  parameters  <7t  are  the  singular  values  of  X^,  u,  and  v,  are  the  corresponding 

left  and  right  singular  vectors,  and  W  is  the  rank  of  X^.   See  [Ref.  38,  Ch.   12]  for 

a  more  detailed  explanation  of  singular  value  decomposition.   Once  a  is  known,  the 

upper  partition  of  (2.2)  can  be  solved  by  simply  carrying  out  the  matrix  multiplication 

described  by  the  left  hand  side  of  (2.5). 

2.      Modified  Yule-Walker  Equation  Methods 

A  well  known  class  of  pole-zero  modeling  techniques  is  based  of  solving 

some  form  of  the  Modified  Yule- Walker  Equations  (MYWE).  These  equations  can  be 

developed  by  multiplying  (2.1)  by  x(n  —  k)  and  taking  the  expectation  of  both  sides. 

This  yields, 

P  Q 

J2  akrxx{n  -k)  =  J2  6trux(n  -  i)  (2.14) 

Jfc=0  t=0 

where  rxx(l)  is  the  autocorrelation  sequence  of  the  system  output  and  rux(l)  is  the 
crosscorrelation  of  the  system  input  and  output.  If  the  original  input  to  (2.1)  is 
assumed  to  be  a  unit  variance  white  noise  sequence  then  the  cross  correlation,  rux(/), 
is  given  by 
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E[u(n)x(n  -  I)]    =    E[u(n)J2h{k)u(n  -  I  -  k)} 

k 

k 

=  H-i). 


(2.15) 


Equation  (2.14)  can  be  then  be  written  as 


p  Q 

J2  akrXx(n  -  k)  =  J2  W*  ~  ")• 

fc=0  t=0 


or  in  matrix  form, 
rxx(0) 

^xx(l) 

rxx(Q) 


rxx{-l) 
rxx(0) 

J'xxi.Q  -]) 

rxx{Q) 


rxx{-P) 
rxx(-P  +  l) 

rxx(Q  -  P) 
'rxx(Q  +  i'-P'j 


rxx{Q  +  P)        rxx(Q  +  P-l) 
.  rxx{Q  +  P+l)        rxx{Q  +  P) 


1 
aP 


rxx{Q  +  1) 
rxx{Q) 

Z?=1bih(i-l) 

bQh(0) 

o 


(2.16) 


(2.17) 


0 

As  in  Prony's  method,  the  solution  for  the  a*  and  6,  coefficients  can  be 
determined  separately.  Taking  a  lower  partition  of  the  last  P  equations  in  (2.17) 
results  in  the  matrix  equation, 

R^a  =  0.  (2.18) 

where  the  theoretical  values  of  the  elements  in  R^  are  replaced  by  estimated  values. 

If  P  +  1  autocorrelation  lags  are  used  in  constructing  R^,  then  (2.18)  can 

be  solved  directly  for  a.  However,  if  additional  reliable  lag  information  is  available,  we 


12 


will  again  desire  to  extend  R^  by  letting  the  index  n  in  (2.16)  run  beyond  Q  +  P  +  1 
resulting  in  additional  equations.  This  leads  to  a  an  overdetermined  set  of  equations, 

R^a  =  e  (2.19) 

which  will  not  in  general  be  satisfied  with  zero  error.  As  before,  application  of  the 
psuedoinverse  results  in  a  least  squares  solution  of  (2.19). 

To  understand  how  the  MYWE  methods  can  be  used  to  match  a  time 
series  to  the  impulse  response  of  a  pole-zero  system  observe  that  (2.16)  describes  the 
relationship  depicted  in  Figure  2.3.  This  operation  can  be  equivalently  expressed  as 

rxx(n)  =  h(n)  *  h{-n).  (2.20) 

If  the  signal  to  be  modeled  is  assumed  to  be  a  pole-zero  system  impulse  response,  then 
for  the  purpose  of  implementing  (2.17),  the  signal  being  modeled  can  be  substituted 
for  h(n)  in  (2.20). 


h(-n) 


B(z 


M 


Figure  2.3:  The  system  relationship  described  by  (2.16). 

The  equations  of  the  upper  partition  of  (2.17)  (the  first  Q  equations)  are 
not  linear  in  the  6  coefficients  and  are  generally  not  solved  directly  from  (2.17).  Once 
the  denominator  coefficients  have  been  determined  from  (2.18)  or  (2.19),  any  of  a 
number  of  techniques  are  available  for  finding  the  desired  transfer  function  numerator 
coefficients.  One  method  already  discussed  is  to  set  up  and  solve  the  upper  partition  of 
Prony's  method  in  (2.5).  Three  other  techniques  are  spectral  factorization,  Durbin's 
method,  and  least  squares  identification,  each  of  which  are  described  below. 
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a.  Spectral  Factorization 

Equation  (2.16)  describes  the  time  domain  relationship 

rxx(l)*a(l)  =  b(I)*h(-l)  (2.21) 

where  a(l)  and  b(l)  represent  the  sequences  of  denominator  and  numerator  transfer 
function  coefficients,  respectively.  Taking  the  2-transform  of  (2.21)  yields 

Sxx(z)A(z)  =  B(z)H(z-x).  (2.22) 

Making  the  substitution 

in  (2.22)  and  rearranging,  results  in 

A(z-l)Sxx(z)A(z)  =  B(z)B(z-1).  (2.24) 

To  utilize  (2.24)  to  find  the  polynomial  coefficients  of  B(z),  which  are  the  elements 
of  the  sequence  b(l),  we  must  perform  spectral  factorization  of  the  sequence  result- 
ing from  the  convolution  of  the  three  sequences  a(/),  a(— /),  and  rxx(l).  A  detailed 
explanation  of  this  technique  can  be  found  in  [Ref.  39]  or  [Ref.  40]. 

b.  Durbin's  Method 

Durbin's  method  [Ref.  41]  makes  use  of  the  property  by  which  a 
process  containing  zeros  can  be  represented  by  an  all-pole  system  if  enough  poles 
are  used.  The  first  step  is  to  filter  the  sequence  to  be  modeled  through  the  inverse 
of  the  previously  determined  all-pole  filter  coefficients  as  shown  in  Figure  2.4.  The 
resulting  residual  sequence,  s(n),  will  nominally  be  an  all-zero  sequence.  A  large 
order  all-pole  model,  Aiarge(z),  can  then  be  fitted  to  s(n)  to  obtain  the  relationship 
illustrated  in  Figure  2.5.  If  the  model  order  for  Aiarge(z)  is  sufficiently  large,  all  of 
the  'information'  in  s(n)  will  be  contained  in  the  coefficients  of  Aiarge(z).   If  an  all 
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x(n) 


s(n) 


Figure  2.4:  Filtering  process  to  generate  the  all-zero  residual  sequence  for 
the  application  of  Durbin's  method. 


u3(n) 


Alargey2) 


s(n) 


Figure  2.5:  The  residual  sequence  approximated  by  a  large  all-pole  model. 

pole  model,  1/B(z),  is  then  constructed  for  the  sequence  of  coefficients  in  AiaTge(z), 
the  relationship  obtained  is 


Alarge{z)  ~ 


1 


B(z) 


(2.25) 


Therefore  the  transfer  function  l/Aiarge(z)  of  Figure  2.5  is  replaced  by 

B(z) 


(2.26) 


Alarge\z) 

which  yields  the  desired  moving  average  (all-zero)  model. 
c.      Shank's  Method 

Consider  the  system  shown  in  Figure  2.6  where  A^(n)  is  the  impulse 
response  of  the  previously  determined  all-pole  portion  of  a  pole-zero  model  (using 
Prony's  method  for  example).  Shank's  method  [Ref.  42]  is  to  satisfy  the  relationship 
of  Figure  2.6  in  the  least  squares  sense.  This  relationship  can  be  described  by  the 
matrix  equation 
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or 


hA(Q)        hA(Q-l) 
hA(Q  +  l)        hA(Q) 

hA(N-l)    hA(N-2) 


hA(0) 

hA(l) 

hA(N-\-Q) 


bo 


>Q 


x(Q) 

x(Q  +  l) 

x(N-l) 


+ 


eB(Q) 

eB(Q  +  l) 

eB(N-l) 


HAb  =  x  +  eB. 


(2.27) 


(2.28) 


Equation  (2.28)  can  be  solved  in  the  manner  of  (2.6)  with  H^  analogous  to  X'^  and  x 
analagous  to  x^.  The  bk  coefficients  can  therefore  be  found  by  using  the  psuedoinverse. 


;(n) 


hA(n) 


Figure  2.6:  System  for  Shank's  method  determination  of  transfer  function 
numerator  coefficients. 


3.      Iterative  Prefiltering 

An  iterative  technique  for  solving  the  direct  modeling  problem  of  Figure  2.1 
called  iterative  prefiltering  has  been  proposed  in  [Ref.  24].  An  effective  application  of 
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the  technique  has  been  reported  in  [Ref.  43].  The  presentation  below  follows  that  of 
[Ref.  40]. 

In  this  method,  the  direct  modeling  problem  error  (see  Figure  2.1a), 

ed(n)  =  x(n)  -  h(n)  (2.29) 

is  expressed  in  the  alternative  form 

€-d{n)    —    x(n)  —  b(n)  *  h,A(n) 

:(n)  *  h,A(n)  *  a(n)  —  b(n)  *  /^(n) 


=     x 


(2.30) 


where  a(n)  and  b(n)  are  the  sequences  of  transfer  function  denominator  and  numer- 
ator coefficients,  respectively,  and  /^(rc)  is  the  impulse  response  of  the  AR  (all-pole) 
portion  of  the  model,  i.e. 

By  then  making  the  equation  error  iterative  (superscripts  represent  the  index  of  iter- 
ation), 

cj+1(n)  =  x(n)  *  tiA{n)  *  ai+l(n)  -  bi+1(n)  *  hA(n),  (2.31) 

the  least  squares  error  solution  for  a1+1(n)  and  bl+l(n)  at  each  iteration  can  be  cal- 
culated using  parameter  estimates  from  the  previous  iteration.  In  matrix  form  (2.31) 
becomes 


*h{P  +  l) 


4(0) 
4(i) 


h\{P) 
hA(P  +  l) 


_xh(N-l)    •••    x\(N-l-P)    hA(N-l) 


hA(P-Q) 
hA(P-Q  +  l) 

hA(N  -1-Q) 


#\P) 
tf\P  +  l) 


1 


a'P+1 


bi+1 
°Q 


(2.32) 
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where  xh(n)  =  x(n)  *  hA(n).  Note  that  (2.32)  can  be  solved  in  exactly  the  same 
manner  as  (2.6). 

a.      Correlation  Domain  Iterative  Prefiltering 

Many  pole-zero  modeling  algorithm's  which  were  originally  conceived 
based  of  the  time  domain  pole-zero  difference  equation,  (2.1),  have  been  reformulated 
in  the  correlation  domain.  Examples  of  this  include  the  correlation  domain  extension 
of  Prony's  method  resulting  in  the  modified  Yule- Walker  equation  methods,  and  an 
instrumental  variable  method  of  least  squares  system  identification  which  Soderstrom 
has  indicated  is  simply  a  correlation  domain  formulation  of  least  squares  system 
identification  [Ref.  44].  In  modeling  trials  conducted  for  this  thesis  both  of  these 
correlation  domain  methods  were  found  to  be  significant  improvements  over  their 
time  domain  counterparts. 

Iterative  prefiltering  can  also  be  extended  into  the  correlation  domain. 
To  see  how  this  is  done  first  note  that  the  direct  formulation  of  the  pole-zero  modeling 
problem  of  Figure  (2.1a)  may  be  reinterpreted  in  the  correlation  domain  by  employing 
the  relationship  of  (2.20)  and  Figure  2.3.  Figure  2.7  illustrates  this  new  direct  pole- 
zero  modeling  interpretation. 

Proceeding  as  for  the  time  domain  iterative  prefiltering  above  we  write 
the  correlation  domain  error  equation, 

ed>COrr(n)    =    rxx(n)  -  x(-n)  *  h(n) 

—    rxx(n)  *  fiA(n)  *  a(n)  —  x(—n)  *  b(n)  *  h,A(n)  (2.33) 

where,  as  before,  a(n)  and  b(n)  are  sequences  of  the  model  coefficients,  /^(n)  is  the 
impulse  response  of  the  AR  (all-pole)  portion  of  the  estimated  model,  and  rxx(n) 
is  found  using  x(n)  as  the  desired  impulse  response  so  that  rxx(n)  =  x(n)  *  x(—n). 
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^xx(tt) 


X 


(-n) 


£<L,corT\?l) 


x(—n)  *  h(n) 


Figure  2.7:  The  direct  pole-zero  modeling  problem  formulation  expressed 
in  the  correlation  domain. 

When  the  error  is  made  iterative,  (2.33)  becomes 


.«+i 


e£r»  =  *M  *  »(-n)  *  h\{n)  *  a*+1(n)  -  «(-n)  *  V+\n)  *  *i(n)  (2.34) 


where  the  superscripts  represent  the  index  of  iteration.  In  matrix  form  (2.34)  becomes 


rh{P) 
rh(P  +  l) 


»4(0) 


x\(P  +  l) 


x\(P  -  Q) 

xi(P-Q  +  l) 


rh(2N-l)    ...    rfaN-l-P)    x\{2N  -  1)    •••    z\(2N  -  1  -  Q) 


<lr(P) 


egSrr(^-l) 


,»'+l 


°o 


(2.35) 


where  r^(n)  =  rxx(n)  *  h\{n)  and  #/i(n)  =  x(—n)  *  h^n)  and  which  can  be  solved  as 
before. 
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4.      Maximum  Likelihood  Techniques 

Maximum  likelihood  estimation  of  parameters  is  the  statistical  standard 
against  which  the  performance  of  other  estimators  is  measured.  This  estimator  makes 
use  of  all  the  useful  statistical  information  available  in  a  given  set  of  data  [Ref.  45,  p. 
73].  Difficulties  arise,  however,  in  the  implementation  of  the  maximum  likelihood  es- 
timator (MLE).  True  maximum  likelihood  estimation  requires  exact  knowledge  of  the 
conditional  probability  density  function  (PDF)  of  the  observed  data  conditioned  on 
the  parameters  to  be  estimated.  This  conditional  PDF  must  then  be  simultaneously 
maximized  with  respect  to  all  parameters  being  estimated.  In  practice,  most  efforts 
to  employ  maximum  likelihood  estimation  make  simplifying  assumptions  about  the 
nature  of  the  input  data  to  derive  a  useful  algorithm.  Such  techniques  are  usually 
called  approximate  maximum  likelihood  methods. 

The  approximate  MLE  chosen  for  this  thesis  is  due  to  Akaike  [Ref.  31]. 
The  brief  developement  of  this  algorithm  provided  below  follows  that  of  Kay  [Ref.  4, 
Ch.  9,10].  Additional  backround  on  maximum  likelihood  estimation  can  be  found  in 
[Ref.  1,  8,  32,  45]  and  references  therein. 

Given  a  sequence  of  independent  random  variable  observations,  x(n),  and 
a  corresponding  set  of  parameters  to  be  estimated,  0*,  the  desired  set  of  estimates  for 
the  0fc's  is  the  one  for  which  the  observed  data  set  is  the  most  likely.  In  terms  of  the 
conditional  probability  density  or  likelihood  function, 

p(x(0),x(l),...,x(7V-l)|^,^2,...,^), 

the  desired  set  of  estimates  is  the  one  which,  for  a  given  set  of  x(n),  is  maximized. 
In  the  case  of  pole-zero  modeling,  an  expression  for  the  observed  data  sequence's 
joint  probability  density  function  conditioned  on  the  model  parameters,  a,  and  bj,  is 
required.  To  obtain  such  an  expression,  it  is  generally  assumed  that  the  observed  data 
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is  Gaussian.  If,  further,  it  is  assumed  that  the  input  to  the  pole  zero  model  is  white 
Gaussian  noise  and  that  the  data  record  length  is  much  longer  than  the  transient 
response  due  to  initial  conditions,  the  conditional  PDF  for  the  observed  data  can  be 
arrived  at  fairly  directly. 

The  joint  PDF  for  the  zero  mean  white  Gaussian  input  sequence,  w(n),  is 
of  the  form 

p(«(0),  u(l),  •  •  • ,  u(N  -  1))  =  n  -L=exp(-^-).  (2.36) 

where  a^  is  the  variance.  Now  the  density  function  for  x(0),  x(l),  •  •  •  ,  x(N  —  1)  con- 
ditioned on  the  0fc's,  can  be  found  from  (2.36)  through  the  standard  linear  transfor- 
mation, 

p(x(0),  *(1),  ■  •  •  ,z(N  -  1))  =  p(u/(0),  11/(1),  •  •  • ,  uf(N  -  1))  \J\ ,  (2.37) 

where  Uf(n)  is  the  inverse  filter  relationship  of  the  original  pole-zero  difference  equa- 
tion, (2.1).  Specifically 

1    p  1    Q 

uf(n)  =  Tz2  aMn  -»)-rE  hu{n  -  k)  (2.38) 

°o  t=0  °o  fc=l 

and  J  is  the  Jacobian  of  the  linear  transformation  uj.  To  simplify  the  transformation 

assume  that  the  pole-zero  filter  in  (2.38)  has  been  expressed  as  its  equivalent  all-pole 

filter, 

Pap 
u(n)  =  ]C  aAP,ix{n  —  i).  (2.39) 

i=0 

The  final  result  of  the  linear  transformation  (2.37)  will  then  be 

p(»(0), a(l), ■  ■  ■  ,x(N  -  1))  =  if  -^L=exp(-£^°  "^f ("  "  k) ).         (2.40) 

n=o  ytncrl  zcru 

When  (2.40)  is  expressed  in  terms  of  the  pole-zero  parameters, 

«ii a2,-'  •  ,ap,bo,bi,--  -,6g, 
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for  the  purpose  of  maximization  the  relationship  is  highly  nonlinear.  Note  that  the 
maximization  of  (2.40)  requires  the  minimization,  over  all  possible  a,'s  and  fcfc's,  of 
the  inverse  filter  error,  uj(n). 

Akaike  employed  the  Newton-Raphson  iterative  method  for  minimizing 
(2.40).  This  method  requires  the  computation  of  Gradient  and  the  Hessian  of  (2.38) 
at  each  iteration  to  generate  the  estimate  updates 


bjfc+i 
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bit 
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dbd&T     ababT 
Using  frequency  domain  arguments,  Akaike  was  able  to  provide  expressions  for  the 

above  partial  derivatives  in  terms  of  Fourier  transforms  which  can  in  turn  be  expressed 

in  terms  of  linear  filtering  operations. 

Note  that  there  are  several  key  assumptions  made  above  which  must  be 

valid  for  this  method  of  approximate  MLE  to  apply: 

1.  The  data  are  real,  Gaussian,  and  zero  mean. 

2.  The  data  record  is  large.  This  is  to  avoid  end  effects  of  assuming  all  data  values 
outside  the  data  record  are  zero  when  filtering  the  data. 

3.  The  poles  and  zeros  are  not  close  to  the  unit  circle.  This  is  to  avoid  long 
transients  due  to  the  initial  conditions  which  are  ignored.  (They  are  assumed 
to  be  known  and  are  set  equal  to  zero.) 

At  first  glance  these  assumptions  would  seem  to  indicate  that  this  method  is  inappro- 
priate to  transient  modeling.  However,  inverse  filter  error  retains  its  meaning  when 
considering  transient  waveforms;  the  ideal  inverse  filtering  result  for  a  transient  signal 
is  a  single  impulse  rather  than  the  minimum  variance  random  sequence  expected  for 
a  stochastic  process.  In  fact,  this  is  exactly  the  appproach  taken  by  Jackson  [Ref.  26, 
pp.  276-278]  in  extending  Judell's  maximum  likelihood  method  [Ref.  33]  to  impulse 
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response  data.  While  this  extension  seems  rather  ad  hoc,  we  will  find  that  such  ap- 
proximate maximum  likelihood  methods  can  be  effective  in  modeling  transient  signals 
as  impulse  responses.  A  key  limitation  imposed  by  dealing  with  deterministic  data 
is  that  reliance  on  inverse  filter  error  excludes  signals  that  must  be  modeled  by  non- 
minimum  phase  systems.  In  contrast,  the  restrictions  of  long  data  record  and  weak 
poles  and  zeros  (not  near  the  unit  circle)  no  longer  apply.  Data  records  end  effects 
and  initial  condition  transient  effects  should  have  no  impact  since  the  assumption  of 
zero  valued  data  outside  the  range  of  data  is  correct  if  the  data  is  chosen  to  end  after 
most  of  the  energy  of  the  transient  is  dissipated. 

C.     IMPLEMENTATION 

All  modeling  algorithms  were  implemented  using  the  interactive  language  PRO- 
MATLAB  from  The  Mathworks,  Inc.  on  SUN  workstations  except  for  the  Akaike 
MLE  algorithm  which  was  implemented  in  FORTRAN  using  a  program  adapted 
from  [Ref.  4,  Ch.  10].  The  FORTRAN  program  was  also  implemented  on  a  SUN 
workstation  with  a  FORTRAN  77  compiler.  All  graphics  were  generated  in  PRO- 
MATLAB. 
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III.  MODELING  PERFORMANCE 

A.     PERFORMANCE  CONSIDERATIONS 

Each  of  the  pole-zero  modeling  techniques  presented  in  Chapter  II  is  effective 
when  modeling  a  signal  that  is  truly  the  impulse  response  of  a  pole-zero  system  with 
no  noise  present  and  with  the  system  order  known.  However,  real  world  transient 
data  rarely  possess  such  characteristics.  Real  world  signals  of  all  types  are  notoriously 
uncooperative  in  fitting  the  signal  models  proposed  to  describe  them.  Reasons  that 
this  may  be  true  for  transient  data  include: 

1.  Inappropriate  selection  of  model  type  or  modeling  algorithm. 

•  Linear  versus  non-linear  models. 

•  Minimum  phase  versus  non-minimum  phase  rational  models. 

2.  The  transient  is  time  shifted  because  of  and  inappropriate  selection  of  the  data 
record  starting  point  due  to  the  presence  of  noise. 

3.  The  assumption  of  impulse  system  excitation  is  a  poor  approximation. 

4.  Incorrect  selection  of  model  order. 

5.  Noise  is  present  in  the  signal. 

The  test  sequences  used  in  this  chapter  are  all  obtained  as  the  impulse  response  of 
linear  pole-zero  systems,  therefore,  the  question  of  linear  versus  non-linear  model  type 
will  not  be  at  issue.  For  the  other  problems,  effective  transient  modeling  requires  both 
selecting  the  appropriate  algorithm  and  understanding  how  to  use  that  algorithm  to 
its  greatest  advantage.    The  next  section  describes  the  test  sequences  used  in  this 
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chapter.    Subsequent  sections  address  how  difficulties  encountered  in  the  pole-zero 
modeling  of  impulse  response  data  can  arise  and  how  they  can  be  overcome. 

The  effects  of  data  selection,  non-minimum  phase  systems,  non-impulse  excita- 
tion, and  incorrect  model  order  on  modeling  will  initially  be  considered  for  signals 
observed  with  no  added  noise  present.  The  effect  of  modeling  a  signal  in  which  addi- 
tive noise  is  present  is  considered  separately.  We  will  see  that  situations  which  modify 
a  linear  pole- zero  system  often  lead  to  another  pole-zero  system.  This  new  system 
usually  has  the  same  number  of  poles  in  the  same  locations  but  with  different  and 
possibly  additional  zeros. 

B.     TEST  SEQUENCES 

The  test  impulse  response  sequences  are  generated  using  pole- zero  models  taken 
from  Kay  [Ref.  4].  The  test  sequence  ARM  A3  uses  one  of  Kay's  models  directly  while 
the  test  sequences  ARMA4  LF,  ARMA3  NM,  and  ARMA4  CL  are  from  Kay  models 
which  have  been  modified  to  enhance  the  illustration  of  certain  points.  The  unit 
impulse  response  and  pole-zero  plot  of  each  test  sequence  model  is  shown  in  Figures 
3.1a-h.  The  model  coefficients  for  these  sequences  are  listed  in  Table  3.1. 


Model 

Model  Coefficients  (a0  =  60  =  1) 

a\ 

<*2 

03 

a\ 

6i 

b2 

ARMA3 

-2.760 

3.809 

-2.654 

0.924 

-0.900 

0.810 

ARMA4  LF 

-3.035 

4.002 

-2.727 

0.778 

-0.200 

0.040 

ARMA3  NM 

-2.760 

3.809 

-2.654 

0.924 

-4.818 

25.000 

ARMA4  CL 

-2.678 

3.700 

-2.5634 

0.917 

-0.200 

0.040 

TABLE  3.1:  Table  of  test  sequence  coefficients. 
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Figure  3.1:  The  test  sequence  ARM  A3,  (a)  Impulse  response  plot  and  (b) 
pole-zero  plot. 
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Figure  3.1:  continued  The  test  sequence  ARMA4  LF.  (c)  Impulse  response 
plot  and  (d)  pole-zero  plot. 
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Figure  3.1:    continued  The  test  sequence  ARM  A3  NM.   (e)  Impulse  re- 
sponse plot  and  (f)  pole-zero  plot. 
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Figure  3.1:  continued  The  test  sequence  ARMA4  CL.  (g)  Impulse  response 
plot  and  (h)  pole-zero  plot. 
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C.     ESTIMATING  NUMERATOR  COEFFICIENTS 

1.      Data  Selection — Time  Shifts  and  Initial  Conditions 

Defining  the  range  of  data  to  be  used  in  modeling  is  an  important  and 
usually  staightforward  exercise.  The  assumption  that  a  signal  represents  the  impulse 
response  of  a  linear  pole-zero  system,  however,  implies  some  very  specific  properties 
about  the  initial  few  points  of  that  signal.  Each  method  of  modeling  the  transfer 
function  numerator  coefficients  in  Chapter  II  reacts  differently  when  the  beginning 
points  of  the  impulse  response  being  modeled  are  degraded.  Because  real  world 
transients  do  not  usually  exhibit  the  instantaneous  rise  time  of  an  ideal  impulse 
response  and  because  noise  is  usually  present,  choosing  the  precise  data  range  for  a 
transient  such  as  that  illustrated  in  Figure  3.2  is  often  a  very  uncertain  task.  Two 
possible  outcomes  when  the  transient  starting  point  is  chosen  incorrectly  are: 

1.  The  starting  point  is  chosen  before  the  signal  begins  so  that  early  data  values 
are  unrelated  (and  presumably  of  lower  amplitude)  to  the  impulse  response  to 
be  matched  (e.g.  these  points  may  consist  of  noise). 

2.  The  starting  point  is  chosen  late  in  which  case  early  values  of  the  impulse 
response  are  lost. 

In  the  first  case,  an  adequate  number  of  additional  model  zeros  can  account 
for  the  delay  in  the  impulse  response.  Assuming  that  the  starting  point  for  the  data  is 
reasonably  close  to  the  true  beginning  of  the  impulse  response,  any  spectral  features 
introduced  by  the  unrelated  early  data  points  will  probably  not  have  enough  energy  to 
significantly  alter  the  spectrum  of  the  impulse  response.  Under  these  circumstances 
all  methods  can  effectively  find  the  poles.  It  is  important  to  note,  however,  that  if 
an  insufficient  number  of  zeros  to  account  for  the  imposed  delay  is  used,  then  some 
of  the  equations  that  are  generated  in  Prony's  method  become  invalid.  When  these 
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Figure  3.2:  An  example  of  a  laboratory  generated  acoustic  transient.  Note 
the  difficulty  in  determining  a  precise  starting  point  for  the  transient. 

equations  axe  solved,  the  invalid  equations  can  drastically  degrade  pole  estimates.  To 
see  this  we  can  apply  Prony's  method  to  a  system  with  the  true  orders  P  =  4  and 
Q  =  2.  Assume  the  signal  is  delayed  by  inserting  three  zeros  at  the  beginning  of  the 
data  so  that  the  original  point  x(0)  is  now  the  fourth  data  point.  If  we  choose  P  =  4 
and  Q  =  3  in  constructing  the  data  matrix  of  (2.2),  the  resulting  set  of  equations  will 
be 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

x(0) 

0 

0 

0 

0 

x(l) 

x(0) 

0 

0 

0 

x(2) 

x(l) 

x(0) 

0 

0 

x(3) 

x(2) 

x(l) 

x(0) 

0 

x(4) 

x(3) 

x(2) 

x(l) 

x(0) 

■  1 

Oi 

b2 

0-2 

ss 

63 

«3 

0 

.  °4  . 

-   *   - 

Note  that  when  the  lower  partition  is  taken  to  find  the  a*  coefficients,  the  first  two 
equations  of  the  lower  partition  are  invalid.   Compounding  the  problem  is  that  the 
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invalid  equations  occur  in  the  high  energy  portion  of  a  transient  signal.  To  overcome 
this  effect  a  numerator  order  of  at  least  five  is  required. 

In  finding  the  numerator  coefficients  of  a  delayed  impulse  response  one 
of  three  outcomes  is  possible  depending  on  the  estimation  technique  chosen.  First, 
any  technique  which  depends  directly  on  the  initial  values  of  the  data  sequence  (for 
example  (2.5))  will  be  ineffective.  Second,  methods  which  rely  on  the  autocorrelation 
of  the  residual  sequence  will  result  in  approximately  the  true  zeros  of  the  system  under 
study.  The  original  time  series  will  not  be  matched  directly.  This  case  is  illustrated 
in  Figure  3.3.  The  resulting  impulse  response  is  an  undelayed  version  of  the  signal. 
Finally,  signal  matching  techniques,  iterative  prefiltering  and  Shank's  method,  result 
in  zeros  not  related  to  the  original  undelayed  model  but  which  provide  the  best  overall 
time  domain  match  of  the  delayed  signal.  This  is  shown  in  Figure  3.4. 

The  case  of  choosing  the  data  record  to  far  to  the  right  and  thus  trun- 
cating the  first  points  of  an  impulse  response  will  again  have  little  effect  on  pole 
estimation.  This  situation  corresponds  the  same  system  with  (non-zero)  initial  con- 
ditions imposed.  Since  initial  conditions  are  acounted  for  in  the  numerator,  the  zeros 
are  significantly  altered.  Here  the  previous  discussion  regarding  finding  the  under- 
lying model  zeros  versus  obtaining  a  good  match  in  the  time  domain  still  applies 
with  one  exception:  the  direct  calculation  of  the  bk  coefficients  from  (2.5)  will  now 
be  effective. 

2.      Non-minimum  Phase  Modeling 

In  [Ref.  46]  it  is  demonstrated  that  the  appropriate  discrete  model  of 
a  sampled  analog  waveform  is  often  represented  by  a  transfer  function  with  zeros 
outside  the  unit  circle.  Many  of  the  techniques  that  are  currently  available  are  pur- 
posefully structured  to  eliminate  uch  non-minimum  phase  models.  In  power  spectrum 
estimation,  a  minimum  phase  system  with  the  same  frequency  response  magnitude 
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as  a  non-minimum  phase  system  results  in  the  same  estimated  spectrum.  Thus  for 
spectrum  estimation  an  equivalent  minimum  phase  system  is  satisfactory.  In  fact, 
the  assumption  underlying  stochastic  modeling  techniques,  namely  white  Gaussian 
noise  input,  guarantees  that  all  transfer  function  combinations  of  minimum  phase 
and  maximum  phase  zeros  are  equivalent.  By  convention,  stochastic  modeling  tech- 
niques always  choose  the  minimum  phase  model  so  that  the  important  statistic  of 
inverse  filter  error  is  available.  In  time  domain  based  applications,  however,  incor- 
rectly choosing  the  model  phase  can  seriously  degrade  system  performance.  Fields 
such  as  seismic  deconvolution,  channel  equilization,  control,  and  matched  filter  de- 
sign, generally  require  identification  of  the  correct  model  phase  [Ref.  47,  48,  49]. 
Also,  we  will  see  in  Chapter  IV  that  effective  modeling  of  real  world  acoustic  signals 
frequently  requires  non-minimum  phase  models. 

A  number  of  modifications  to  the  basic  stochastic  model  have  been  in- 
troduced to  allow  selection  of  the  model  with  the  correct  phase.  These  techniques 
generally  involve  changing  the  Gaussian  nature  of  the  input  noise  [Ref.  50]  and  often 
employ  higher  than  second  order  moments  or  cumulants  [Ref.  29,  30].  However,  when 
a  model's  impulse  response  has  effectively  matched  a  signal  in  the  time  domain,  the 
resulting  phase  is  immediately  known  to  be  the  correct.  Allowing  for  the  possibility 
of  non-minimum  phase  models  places  significant  limitations  on  the  modeling  methods 
which  may  be  used.  Techniques  which  rely  on  inverse  filtering  (maximum  liklihood 
methods)  are  not  applicable  since  the  inverse  of  a  non-minimum  phase  system  is  unsta- 
ble. Also,  techniques  which  use  correlation  information  to  calculate  the  bk  coefficients 
(spectral  factorization  and  Durbin's  method)  will  give  poor  results  since  correlation 
data  does  not  preserve  phase  information.  Figure  3.5  shows  the  difficulty  encountered 
when  Durbin's  method  attempts  to  model  a  non-minimum  phase  system.  The  best 
Durbin's  method  can  do  is  produce  the  spectrally  equivalent  minimum  phase  version 
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of  a  maximum  phase  system  since  the  autoregressive  modeling  techniques  on  which 
it  relys  can  only  produce  zeros  within  the  unit  circle.  In  contrast,  Figure  3.6  shows 
that  Shank's  method  is  able  to  find  the  correct  model. 
3.      Numerator  Modeling  Summary 

Table  3.2  provides  a  brief  summary  of  the  modeling  properties  of  the  nu- 
merator coefficient  modeling  techniques  considered  in  this  thesis.  Since  the  goal  in 
Chapter  IV  is  to  perform  time  domain  modeling  of  the  acoustic  transients  being 
considered,  those  methods  which  provide  the  best  time  domain  match  between  the 
original  signal  and  the  model  impulse  response  are  preferred. 


Method 

Equation 

Non-minimum 
phase  capable? 

Data  Selection 

Shifted  Right 
(Delayed) 

Shifted  Left 
(Truncated) 

Prony,  upper 
partition 

2.5 

Yes 

Not 
Usable 

Time  Series 
Match 

Spectral 
Factorization 

2.24 

No 

Underlying 
Model 

Underlying 
Model 

Durbin's 
Method 

2.25 
2.26 

No 

Underlying 
Model 

Underlying 
Model 

Shank's 
Method 

2.28 

Yes 

Time  Series 
Match 

Time  Series 
Match 

Iterative 
Prefiltering 

2.32 

Yes 

Time  Series 
Match 

Time  Series 
Match 

Akaike 
MLE 

2.41 

No 

Underlying 
Model 

Time  Series 
Match 

TABLE  3.2:    Summary  of  the  capabilities  and  limitations  of  numerator 
modeling  methods. 


D.     NON-IMPULSE  EXCITATION 

In  any  real  world  system  the  assumption  of  a  unit  impulse  input  is  approxi- 
mate. If  the  duration  of  the  excitation  waveform  is  small  relative  to  the  period  of  the 
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Figure  3.3:  The  right  shifted  sequence  ARMA3,  poles  modeled  using 
LSMYWE  and  zeros  modeled  using  Durbin's  method,  (a)  Time  signal 
plot  and  (b)  pole-zero  plot.  Durbin's  method  does  not  account  for  the 
time  delay  but  instead  finds  the  underlying  system's  true  zeros. 
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Figure  3.4:  The  right  shifted  sequence  ARMA3,  poles  modeled  using 
LSMYWE  and  zeros  modeled  using  Shank's  method,  (a)  Time  signal 
plot  and  (b)  pole-zero  plot.  The  least  squares  method  does  not  find  the 
underlying  system's  true  zeros  but  rather  achieves  the  best  overall  time 
domain  match. 
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Figure  3.5:  The  non-minimum  phase  sequence  ARMA3  NM,  poles  mod- 
eled using  LSMYWE  and  zeros  modeled  using  Durbin's  method,  (a)  Time 
signal  plot  and  (b)  pole-zero  plot.  Durbin's  method  cannot  model  zeros 
outside  the  unit  circle. 


37 


original  signal 

model  impulse  response 


(a) 


x  true  poles 

o  true  zeros 

1.5 


150    200 


+  estimated  poles 
*  estimated  zeros 


(b) 


z-plane 

Figure  3.6:  The  non-minimum  phase  sequence  ARMA3  NM,  poles  mod- 
eled using  LSMYWE  and  zeros  modeled  using  Shank's  method,  (a)  Time 
signal  plot  and  (b)  pole-zero  plot.  The  least  squares  method  is  effective 
at  modeling  zeros  outside  the  unit  circle. 
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lowest  oscillation  frequencies  present  then  the  assumption  is  justified.  However,  it  is 
reasonable  to  expect  this  criterion  will  often  not  be  satisfied.  Non-impulse  excitations 
which  may  be  encountered  include: 

1.  A  long  duration  baseband  type  pulse. 

2.  An  uncorrected  random  train  of  impulses.  This  model  is  often  used  to  account 
for  reverberation  (echoes)  in  siesmic  deconvolution  [Ref.  51]. 

3.  A  frequency  swept  input  that  sweeps  through  the  natural  frequencies  of  a  sys- 
tem. This  model  is  usually  considered  in  conjunction  with  the  starting  and 
stopping  of  rotating  machinery. 

If  the  system  input  were  known,  the  modeling  problem  could  be  formulated  as  a  sys- 
tem identification  problem.  When  no  information  about  the  system  input  is  available, 
other  means  must  be  found  to  deal  with  this  problem. 
1.     Baseband  Pulse  Excitation 

The  effect  of  modeling  a  transient  signal  from  a  linear  system  in  which  the 
input  is  a  long  duration,  baseband-type  pulse  can  best  be  understood  as  filtering  by  a 
finite  impulse  response  (FIR,  all-zero)  filter  as  illustrated  in  Figure  3.7.  The  spectral 
properties  of  the  original  time  series  are  windowed  by  the  frequency  response  of  the 
FIR  filter  coefficients.  For  this  type  of  pulse,  the  effect  is  that  of  low  pass  filtering. 
Thus  high  frequencies  are  attenuated  relative  to  low  frequencies.  If  no  frequency 
component  exists  below  the  FIR  filter's  cutoff  frequency  then  the  original  spectrum 
is  altered  according  to  the  side  lobe  structure  of  the  FIR  filter. 

The  new  model  that  results  can  be  viewed  as  a  system  with  the  original 
model  poles  but  with  new  numerator  polynomial  coefficients  that  are  the  result  of 
convolving  the  FIR  filter  coefficients  with  the  original  numerator  polynomial  coeffi- 
cients. This  results  in  a  higher  order  polynomial,  hence  more  zeros  than  were  present 
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Figure  3.7:  One  way  to  view  a  linear  system  excited  by  a  baseband  pulse. 

in  the  original  model  will  be  required.  Observe  how  the  original  impulse  response  of 
Figure  3.7a  is  altered  to  Figure  3.7b  when  a  nine  element  triangular  excitation  pulse 
is  used.  The  dotted  lines  in  Figures  3.7b,c,d  indicate  the  modeling  results  obtained 
when  none,  two,  and  four  extra  zeros,  respectively,  are  used  in  the  estimated  pole-zero 
model.  In  this  case  four  extra  zeros  prove  sufficient  to  account  for  the  input  pulse. 
The  corresponding  model  spectrum,  Figure  3.7e  illustrates  the  attenuation  caused  by 
the  baseband  excitation.  The  pole-zero  plot  in  Figure  3.7f  shows  that  non-minimum 
phase  zeros  were  required  to  achieve  an  effective  time  domain  match. 

2.  Random  Impulse  Train  Excitation 

If  the  input  to  a  linear  system  is  an  uncorrelated  random  train  of  impulses 
then,  although  the  time  series  may  be  significantly  different  from  the  original  impulse 
response,  the  autocorrelation  function  of  the  signal  is  theoretically  unaltered  except 
for  a  scaling  factor.  This  is  because  the  autocorrelation  of  an  uncorrelated  impulse 
train  is  a  scaled  unit  impulse.  Thus  modeling  methods  which  rely  on  correlation 
information  should  be  effective.  However,  over  a  finite  time  interval  it  is  unlikely  that 
a  random  sequence  will  be  truly  uncorrelated.  As  the  impulse  train  becomes  correlated 
the  situation  will  be  equivalent  to  the  baseband  pulse  case  described  above. 

3.  Frequency  Swept  Excitation 

The  output  of  a  system  excited  with  a  frequency  swept  signal  depends  on 
the  rate  at  which  the  sweep  occurs.  A  slow  sweep  will  result  in  a  series  of  transient 
events,  each  at  a  specific  resonant  frequency  of  the  system.   These  events  can  each 
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be  modeled  separately.  When  the  sweep  rate  is  rapid,  all  natural  modes  will  appear 
much  as  if  the  input  were  an  impulse  except  that  the  phase  relationships  of  the  various 
components  may  be  changed. 

E.     MODEL  ORDER  SELECTION 

In  studies  of  rational  modeling  the  issue  which  continues  to  be  the  most  con- 
founding is  that  of  model  order  selection.  The  proposed  methods  which  have  a  sound 
theoretical  basis  (e.g.  [Ref.  52,  53,  54])  are  very  difficult  to  actually  implement. 
These  methods  are  invariably  related  to  maximum  likelihood  concepts  and  therefore 
rely  heavily  on  inverse  filter  statistics.  This  implies  that  for  model  order  evaluation 
the  inverse  filter  error  must  be  calculated  over  all  possible  model  orders.  Then  the 
model  order  and  inverse  filter  error  which  minimize  some  function  of  the  two  is  se- 
lected. The  case  of  non-minimum  phase  systems  is  even  more  intractable  since  the 
inverse  filter  of  such  systems  is  unstable. 

A  more  attractive  but  less  understood  method  is  to  initially  overdetermine  a 
system  and  then  allow  the  modeling  algorithm  to  indicate  the  correct  model  order. 
A  method  proposed  by  Cadzow  [Ref.  39]  uses  singular  value  decomposition  to  aid  in 
model  order  selection  in  the  denominator  of  all-pole  and  pole-zero  models.  Kumeresan 
and  Tufts  [Ref.  21]  have  shown  that  when  Prony's  method  is  applied  to  exponentially 
damped  sinusoids  reversed  in  time,  valid  poles  occur  outside  the  unit  circle  and  excess 
poles  occur  inside  the  unit  circle.  In  both  of  these  methods  the  denominator  order 
is  initially  overdetermined  and  the  modeling  method  then  provides  the  correct  order. 
For  this  thesis  both  of  these  methods  were  applied  to  a  number  of  transients.  They 
were  effective  when  applied  to  the  noiseless  impulse  responses  of  true  pole- zero  systems 
but  were  not  robust  in  the  presence  of  noise  or  when  many  narrowband  components 
are  present.  There  is  no  similar  guidance  for  determining  the  numerator  model  order. 
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Figure  3.8:  The  sequence  ARMA4  LF  with  (a)  unit  impulse  excitation 
and  (b)  triangular  pulse  input  modeled  with  correct  model  orders  P  —  4 
and  Q  =  2  using  Prony's  method  and  Shank's  method. 
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Figure  3.8:  continued  The  sequence  ARMA4  LF  with  triangular  pulse 
input  (c)  modeled  with  Prony's  method  and  least  squares  identification 
with  orders  P  =  4  and  Q  =  4  (d)  modeled  with  Prony's  method  and  Shank's 
method  orders  P  =  4  and  Q  =  6. 
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Figure  3.8:  continued  The  sequence  ARMA4  LF  with  triangular  input  (e) 
the  original  impulse  response  and  triangular  output  model,  order  (4,6), 
spectrums  and  (f)  the  corresponding  pole-zero  plot.  Note  how  the  base- 
band pulse  has  attenuated  the  higher  frequencies. 
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Perhaps  most  desirable  would  be  a  modeling  technique  that,  when  overdeter- 
mined,  caused  excess  poles  and  zeros  to  either  cancel  or  migrate  well  away  from  the 
unit  circle  (all  the  way  to  the  origin  ideally).  Tummala  [Ref.  23]  has  demonstrated 
some  success  with  this  concept  using  an  iterative  algorithm  to  solve  the  least  squares 
identification  problem.  Observing  the  degree  to  which  the  various  transient  modeling 
techniques  of  this  thesis  exhibit  this  behavior  is  the  approach  taken  in  the  following 
comparison.  This  behavior  was  observed  by  modeling  the  ARMA3  test  sequence  with 
different  combinations  of  numerator  and  denominator  order.  Table  3.3  summarizes 
the  results  obtained.  A  'Y'  in  Table  3.3  indicates  a  modeling  technique  achived  an  ex- 
act time  domain  match  between  the  model  impulse  response  the  modeled  signal.  The 
most  notable  negative  result  is  that  both  Prony's  method  and  LSMYWE  were  ineffec- 
tive when  the  numerator  and  denominator  were  overdetermined.  Also,  excess  poles 
without  enough  zeros  causes  problems  for  correlation  domain  iterative  prefiltering. 
Figure  3.9  is  an  example  of  the  results  obtained  using  LSMYWE  when  excess  poles 
and  zeros  are  present.  However,  when  these  results  were  used  to  initialize  iterative 
prefiltering  the  excess  poles  and  zeros  were  handled  effectively. 

The  above  results  and  experience  gained  during  extensive  modeling  trials  lead  to 
the  following  recommended  strategy  when  modeling  complex  transients  about  which 
very  little  is  known: 

1.  Be  conservative  in  estimating  the  denominator  order.  Signals  composed  of  nar- 
rowband components  will  generally  be  dominated  by  the  few  components  high- 
est in  energy.  Even  if  excess  zeros  are  required  to  deal  with  a  noisy  signal  (see 
the  next  section)  a  small  number  of  excess  zeros  will  usually  suffice  when  us- 
ing LSMYWE.  Pole  detection  techniques  such  as  those  of  Cadzow[Ref.  5],  and 
Kumeresan,  and  Tufts  [Ref.  21]  may  be  helpful. 
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Method 

Correctly  Modeled?  (Y/N) 

P=10 
Q=2 

P=4 
Q=10 

P=10 
Q=4 

P=6 
Q=10 

Prony's 
Method 

Y 

Y 

N 

N 

LSMYWE 

Y 

Y 

N 

N 

Iterative 
Prefiltering 

Y 

Y 

Y 

Y 

Corr  Domain 
Iterative  Pref 

N 

Y 

Y 

Y 

Akaike 
MLE 

Y 

Y 

Y 

Y 

TABLE  3.3:  The  effectiveness  of  thesis  modeling  methods  on  ARMA3 
with  varying  overdetermination  of  model  order.  'Y'  indicates  an  exact 
time  domain  match  was  achieved  and  *N'  indicates  a  poor  time  domain 
match  resulted. 

2.  Be  expansive  with  estimates  of  numerator  order.  (Iterative  prefiltering  may  be 
an  exception  to  this.  See  Chapter  IV.)  Additional  zeros  can  help  compensate 
for  mistakes  made  in  data  selection  and  effects  of  non-impulse  excitation.  Un- 
needed  zeros  will  usually  migrate  away  from  the  unit  circle.  Note  that  when 
using  Prony's  method  or  LSMYWE,  one  numerator  order  can  be  used  when 
calculating  the  denominator  coefficients  and  another  numerator  order  can  be 
chosen  when  finding  the  numerator  coefficients.  Using  Shank's  method  with 
Prony's  method  or  LSMYWE  provides  the  most  flexibility  when  a  model  is 
attempting  to  account  for  erroneous  assumptions  (e.g.  data  selection)  which 
may  have  been  made  by  the  user.  The  above  recomendations  are  primarily 
intended  for  Prony's  method  and  LSMYWE.  However,  using  any  iterative  tech- 
nique requires  reasonable  initial  estimates  which  are  usually  arrived  at  using 
linear  estimation  techniques. 
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Figure  3.9:  Pole-zero  plots  for  the  sequence  ARMA3  modeled  with  overde- 
termined  model  orders  of  P=10  and  Q=10.  (a)  lsmywe  and  Durbin's 
method  results  in  an  unstable  system,  (b)  iterative  prefiltering  provides 
an  effective  estimation  of  the  correct  poles  and  zeros  and  an  excellent  time 
domain  match  (not  shown). 
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F.     NOISE  PERFORMANCE 

Kay  [Ref.  55]  has  shown  that  for  an  all-pole  process,  additive  white  noise  has 
the  effect  of  introducing  zeros  which  migrate  from  the  origin  to  the  model  poles  as 
the  signal-to-noise  ratio  is  decreased.  Alternatively,  if  zeros  are  not  introduced  into 
the  model,  poles  move  toward  the  origin  as  the  signal-to-noise  ratio  is  reduced  [Ref. 
56].  In  either  case  the  overall  effect  is  a  loss  of  spectral  resolution.  This  result  extends 
directly  to  pole-zero  processes  in  white  noise.  Therefore  one  important  feature  of  any 
modeling  technique  is  it's  ability  to  resolve  spectral  components  in  the  presence  of 
noise. 

1.      Rational  Modeling  of  Noisy  Data 

Several  authors  have  recently  addressed  the  difficulties  associated  with 
modeling  a  time  series  in  which  additive  white  noise  is  present.  Kay  [Ref.  57]  has 
noted  that  when  the  variance  of  the  white  noise  can  be  estimated,  its  effect  can 
be  removed  from  the  zeroth  autocorrelation  lag  to  improve  the  resolution  of  pole 
frequencies  in  correlation  based  autoregressive  modeling.  Alternatively,  the  zeroth 
autocorrelation  lag  may  be  avoided  by  using  (2.18),  sometimes  called  the  high  order 
Yule- Walker  Equations,  for  Q  >  P  or  by  eliminating  the  rows  containing  the  zeroth 
lag  in  certain  least  squares  formulations  [Ref.  44].  If  the  resulting  high  order  Yule- 
Walker  equations  yield  a  matrix  that  is  well  conditioned,  this  procedure  is  equivalent 
to  that  of  Kay  above  [Ref.  58].  Even  with  good  correlation  estimates,  however,  the 
matrix  R^  is  not  guaranteed  to  be  positive  definite  (i.e.,  invertible).  Consequently, 
R4  is  more  likely  to  be  a  poorly  conditioned  matrix.  This  is  not  a  problem  in  the 
least  squares  formulation,  (2.19),  since  the  very  act  of  choosing  such  a  formulation 
implies  we  expect  a  solution  that  is  approximate. 

A  number  of  authors  have  noted  that  overdetermining  the  number  of  model 
poles  can  have  a  profound  effect  when  modeling  noisy  signals  [Ref.  55,  56,  21,  36,  39]. 
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The  explanation  for  this  is  that  the  extra  poles  are  able  to  model  the  flat  noise 
spectrum  of  the  signal,  preventing  that  portion  of  the  signal  from  interfering  with 
the  poles  needed  for  narrowband  modeling.  Finally,  singular  value  decomposition  has 
been  used  to  aid  resolution  by  overcoming  the  conditioning  problems  that  can  result 
from  pole  overdetermination  and  high  order  Yule- Walker  equations  [Ref.  5,  21,  36, 
59,  60]. 

Little  work  has  appeared  which  analyzes  or  demonstrates  the  effectiveness 
of  iterative  prefiltering  in  the  presence  of  additive  white  noise.  Stoica  and  Soder- 
strom  [Ref.  61]  have  shown  that  iterative  prefiltering  will  perform  well  for  very  long 
stochastic  data  sequences  in  the  presence  of  additive  white  noise  if  reasonable  initial 
estimates  are  available  to  start  the  iterations.  They  also  show  that  for  arbitrary  initial 
estimates  iterative  prefiltering  is  not  guaranteed  to  converge.  Tufts  and  Kumeresan 
[Ref.  62]  have  demonstrated  superior  performance  for  approximate  MLE  over  linear 
prediction  methods  when  modeling  sinusoids  in  white  noise.  Using  the  approximate 
MLE  method  of  Box  and  Jenkins,  Cadzow  found  that  such  a  method  performed  com- 
parably to  the  overdetermined  high  order  Yule- Walker  equations  for  the  sum  of  two 
second  order  all-pole  processes  in  white  noise  but  with  higher  variance. 

The  above  work  suggests  numerous  possibilities  in  dealing  with  noisy  im- 
pulse response  data.  The  effectiveness  of  these  techniques  is  evaluated  empirically  in 
the  following  section  by  modeling  noisy  test  sequences. 
2.     Discussion  and  Examples 

The  examples  used  to  illustrate  modeling  performance  are  summarized 
in  Table  3.4.  All  noise  modeling  examples  were  conducted  using  the  test  sequence 
ARMA4  CL  (see  Table  3.1).  Figure  3.10  shows  that  Prony's  method  has  difficulty 
even  when  the  SNR  is  as  high  as  30  dB  although  the  excess  poles  introduced  for 
Figure  3.10b  dramatically  increase  modeling  effectiveness.  The  moderate  impact  of 
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insuring  that  Q  >  P  is  illustrated  for  LSMYWE  in  Figure  3.11  in  that  the  higher 
frequency  pole  moves  out  toward  the  unit  circle  when  the  modeling  order  is  altered 
from  (4,2)  to  (4,4).  As  with  Prony's  method,  overestimating  the  number  of  poles 
dramatically  improves  modeling  in  Figures  3.12  and  3.13.  In  Figure  3.14,  the  two 
previous  noise  modeling  examples  are  again  modeled  with  the  smallest  singular  value 
set  to  zero  to  yield  a  data  matrix  of  rank  P.  The  technique  is  highly  effective  in 
both  cases.  The  effectiveness  of  singular  value  decomposition  when  overdetermining 
the  number  of  poles  is  illustrated  by  the  examples  in  Figure  3.15.  This  is  identical 
to  the  result  described  in  [Ref.  21].  Discarding  the  excess  singular  values  aids  both 
in  resolving  the  narrowband  components  present  and  in  reducing  the  variance  of  the 
excess  poles.  Figures  3.16  and  3.17  show  the  hazards  that  are  possible  when  using 
singular  value  decomposition.  Any  number  of  singular  values  discarded  other  than 
the  those  necessary  to  achieve  a  data  matrix  rank  equal  to  the  transfer  function 
denominator  order  has  undesirable  side  effects.  Discarding  too  few  singular  values 
such  as  in  the  Figure  3.16  examples  will  resolve  the  desired  narrowband  components 
but  also  lead  to  excess  poles  near  or  even  beyond  the  unit  circle.  This  result  will 
at  best  increase  the  likelihood  of  spurious  narrowband  components  and  can,  in  the 
worst  case,  lead  to  an  unstable  model.  If  too  many  singular  values  are  discarded 
as  in  Figure  3.17b,  the  accuracy  of  narrowband  component  estimates  will  be  badly 
degraded.  Therefore,  singular  value  decomposition,  while  potentially  very  useful  in 
modeling  transient  signals  in  noise,  must  be  used  with  caution  when  the  true  model 
order  is  unknown  (which  is  basically  always). 

The  iterative  prefiltering,  correlation  domain  iterative  prefiltering,  and 
Akaike  MLE  algorithms  were  initialized  with  several  of  the  preceeding  examples  and 
the  results  are  shown  in  and  in  Figures  3.18  and  3.19.  Iterative  prefiltering  dramati- 
cally improved  the  narrowband  modeling  of  each  example  with  which  it  was  initialized. 
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Figure  3.18a  indicates  that  excess  poles  may  cause  problems  by  being  close  to  the  unit 
circle.  However,  the  ability  of  iterative  preflltering  to  cancel  excess  zeros  with  excess 
poles  demonstrated  in  this  chapter  may  mitigate  this  effect.  Correlation  domain  iter- 
ative preflltering  also  improved  the  resolution  of  each  noise  example  on  which  it  was 
tried  although  its  convergence  was  slower  than  and  the  bias  of  its  results  were  higher 
than  standard  iterative  preflltering.  Also,  correlation  domain  iterative  preflltering 
showed  a  alarming  tendency  to  place  excess  poles  near  the  true  poles  on  noisy  sig- 
nals. We  experienced  a  great  deal  of  difficulty  in  trying  to  model  noisy  signals  with 
the  Akaike  MLE  algorithm.  The  primary  problem  encountered  was  early  algorithm 
termination  due  to  an  iteration  which  produced  a  non-minimum  phase  zero.  Akaike 
MLE  was  able  to  improve  the  quality  of  a  'good'  Prony  estimate  in  Figure  3.19b. 
Finally,  Figure  3.20  shows  a  pole-zero  modeling  result  corresponding  to  the  example 
in  Figure  3.17a.  Figure  3.20a  shows  the  noisy  signal  used  to  during  modeling  and 
Figure  3.20b  compares  the  resulting  model  to  the  original  noiseless  impulse  response. 

G.     MODELING  PERFORMANCE  SUMMARY 

In  modeling,  the  more  that  is  known  about  the  signal  to  be  modeled,  the  better 
the  choice  of  model  structure  and  modeling  algorithm  that  can  be  made.  The  focus 
of  this  thesis  is  the  modeling  of  real  world  transient  signals  about  which  very  little  is 
known  and  whose  characteristics  can  vary  widely.  The  type  of  model  is  also  set  by  the 
basic  goal  of  this  thesis,  that  is,  a  rational  linear  model.  Under  such  circumstances, 
a  sensible  criterion  in  choosing  a  modeling  algorithm  is  robustness.  For  numerator 
modeling,  Shank's  method  is  particularly  robust  in  the  sense  that  it  provide  the 
numerator  coefficients  which  give  the  minimum  norm,  least  squares  fit  based  on  the 
previously  determined  denominator  coefficients.  Shank's  method  is  insensitive  to 
time  shifted  data  records  and  uncertainties  in  model  order.  Durbin's  method  is  also 
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effective  but  suffers  from  the  limitations  that  only  minimum  phase  models  are  possible 
and  the  resulting  models  cannot  account  for  time  shifts.  Since  the  Akaike  MLE 
algorithm  requires  a  minimum  phase  initial  estimate,  Durbin's  method  can  be  used 
provide  an  such  an  estimate.  Equation  (2.5)  and  spectral  factorization  are  extremely 
sensitive  to  degraded  signals  (e.g.  noisy  signals)  and  so  will  not  be  considered  further. 
Denominator  modeling  with  Prony's  method  is  extremely  sensitive  to  even  small 
amounts  of  noise.  This  is  only  slightly  less  true  of  Akaike  MLE.  Also,  Akaike  MLE 
cannot  continue  to  iterate  if  a  non-minimum  phase  model  estimate  is  encountered.  In 
contrast,  LSMYWE  and  iterative  prefiltering  are  quite  robust  in  noise  and  iterative 
prefiltering  is  particularly  robust  to  model  order  overdetermination.  Our  conclusion 
is  that  for  modeling  the  real  world  acoustic  transients  of  the  next  chapter,  LSMYWE 
with  Shank's  method  and  iterative  prefiltering  are  likely  to  perform  the  best  on  the 
acoustic  transients  considered  in  the  next  chapter. 
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Method(order  P,Q) 
(initialization) 

SNR 
in  dB 

Bias                   Variance  of  Pole 
(true  minus  model)            Estimates 

Figure 

True  Poles 

.6247  ±  .7509; 
.7144  ±. 6711; 

Prony  (4,2) 

30 

-.0598  ±.0330; 
poor 

<  10~4 
poor 

3.10a 

Prony  (6,2) 

30 

.0283  ±. 0186; 
.0176 +. 0074; 

<  10"4 

<  10"4 

3.10b 

LSMYWE  (4,2) 

15 

. 1788  ±. 1872; 
.0346 +  .0180; 

<  10"4 
.0012 +.0014; 

3.11a 

LSMYWE  (4,4) 

15 

.1475  ±.1035; 
.0386  =F  .0045; 

<  10"4 
.0022  +.0019; 

3.11b 

LSMYWE  (6,2) 

15 

.0053  ±  .0035; 
.0027  ±  .0014; 

<  10"4 

<  10"4 

3.12a 

LSMYWE  (8,8) 

10 

.0040  ±  .0054; 
.0002  ±  .0037; 

<  10"4 

<  10"4 

3.12b 

LSMYWE  (8,8) 

5 

.03 16  ±. 0497; 
.0147  ±. 0212; 

.0002 +.0001; 
.0008  +  .0010; 

3.13a 

LSMYWE  (12,12) 

5 

.0115  ±. 0077; 
.0017  ±. 0063; 

<  10"4 

<  10"4 

3.13b 

Prony  (4,2)  SVD 

30 

-. 0010  ±.0040; 
-.0016  ±.0014; 

.0001  +  .0002; 
<  10~4  +  .003; 

3.14a 

LSMYWE  (4,4)  SVD 

15 

.0017  T  0005; 
-.0007  ±.0005; 

<  10"4 

<  10"4 

3.14b 

Prony  (8,2) 

20 

-. 0030  ±. 0387; 
.0275  ±. 0188; 

.0001  +  .0001; 
.0001  +  .0001; 

3.15a 

Prony  (8,2)  SVD 

20 

-.0105  ±. 0043; 
-. 0059  ±. 0077; 

<  10"4 

<  10~4 

3.15b 

LSMYWE  (8,8)  SVD 

5 

-.0074  ±.0100; 
-.0070  ±.0086; 

.0006  +  .0008; 
.0003  +  .0009; 

3.16a 

LSMYWE  (8,8)  SVD 

5 

-.0012  ±.0069; 
-. 0062  ±. 0064; 

.0004  +  .0003; 
.0004  +  .0008; 

3.16b 

LSMYWE  (8,8)  SVD 

5 

-.0037  ±.0024; 
-.  0049  ±.  0022; 

.0001  +  .0001; 
.0002  +  .0005; 

3.17a 

LSMYWE  (8,8)  SVD 

5 

-.0342  ±.0526; 
-. 0507  ±. 0730; 

.0001  +  .0002; 
<  10"4 

3.17b 

Iterative 
Prefiltering  (4,4) 

15 

-.0001  ±. 0001; 
.0002  ±.0001; 

<  10"4 

<  10~4 

3.18a 

Iterative 
Prefiltering  (8,8) 

5 

-.0001  ±.0011; 
-.0006  ±. 0002; 

<  10"4 

<  10~4 

3.18b 

Corr  Domain  (8,8) 
Iter  Pref 

5 

-.0054  ±.0030; 
.0002  ±  .0017; 

<  10"4 

<  10"4 

3.19a 

Akaike  MLE  (6,2) 

30 

-.0008 +.0009; 
-.0004  t. 0010; 

<  10"4 

<  10"4 

3.19b 

TABLE  3.4:  Examples  of  pole-zero  modeling  performance  of  the  impulse 
response  test  sequence  ARMA4  CL  in  added  noise. 
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Figure  3.10:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  ef- 
fectiveness of  overdetermining  the  number  of  model  poles  using  Prony's 
method.  Test  sequence  is  ARMA4  CL  with  30  dB  of  noise,  (a)  Using  the 
correct  model  order,  (4,2)  and  (b)  Using  two  excess  poles,  model  order 
(6,2). 
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Figure  3.11:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  mod- 
erate benefits  of  zeroth  lag  correlation  compensation  by  choosing  Q  —  P 
when  using  LSMYWE.  Test  sequence  is  ARMA4  CL  with  30  dB  of  noise, 
(a)  Using  the  correct  model  order,  (4,2)  and  (b)  model  order  (4,4). 
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Figure  3.12:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the 
dramatic  benefits  of  using  excess  poles  with  LSMYWE.  Test  sequence 
is  ARMA4c.  (a)  Model  order  (6,2)  with  15  dB  of  added  noise  and  (b) 
model  order  (8,8)  with  lOdB  of  added  noise. 
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Figure  3.13:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the 
dramatic  benefits  of  using  excess  poles  with  LSMYWE.  Test  sequence 
is  ARMA4  CL.  (a)  Model  order  (8,8)  with  5  dB  of  added  noise  and  (b) 
model  order  (12,12)  with  5  dB  of  added  noise. 
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Figure  3.14:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  dra- 
matic benefits  of  setting  the  smallest  singular  value  of  the  data  matrix  to 
zero  (i.e.  adjust  data  matrix  rank  to  P)  for  Prony's  method  and  LSMYWE. 
Test  sequence  is  ARMA4  CL.  (a)  Prony's  method  for  model  order  (4,2) 
with  30  dB  of  added  noise  and  (b)  LSMYWE  for  model  order  (4,4)  with 
15  dB  of  added  noise. 


58 


x  true  poles 

.  estimated  poles 


0.5  - 


-0.5  - 


(a) 


z-plane 


x  true  poles 

.  estimated  poles 


0.5  - 


-0.5  - 


(b) 


z-plane 


Figure  3.15:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  effect 
of  adjusting  the  data  matrix  rank  when  using  excess  poles  for  Prony's 
method.  Test  sequence  is  ARMA4c.  (a)  Model  order  (8,2)  with  20  dB 
of  added  noise  and  no  rank  adjustment  and  (b)  model  order  (8,2)  with 
20  dB  of  added  noise  with  rank  adjusted  to  Ptrue  =  4  using  singular  value 
decomposition. 
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Figure  3.16:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  effect 
of  adjusting  the  data  matrix  rank  when  using  excess  poles  for  LSMYWE. 
Test  sequence  is  ARMA4  CL.  (a)  Model  order  (8,8)  with  5  dB  of  added 
noise  with  rank  adjusted  to  Ptme  +  4  =  8  and  (b)  model  order  (8,8)  with  5 
dB  of  added  noise  with  rank  adjusted  to  PtTue  +  2  =  6. 
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Figure  3.17:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  effect 
of  adjusting  the  data  matrix  rank  when  using  excess  poles  for  LSMYWE. 
Test  sequence  is  ARMA4  CL.  (a)  Model  order  (8,8)  with  5  dB  of  added 
noise  with  rank  adjusted  to  Ptrxit  =  4  and  (b)  model  order  (8,8)  with  5  dB 
of  added  noise  with  rank  adjusted  to  PtTue  —  1  =  3. 
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Figure  3.18:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  abil- 
ity of  iterative  prefiltering  to  improve  the  resolution  of  an  LSMYWE  es- 
timate. In  both  cases  the  iterative  prefiltering  algorithm  was  initialized 
using  LSMYWE  and  Durbin's  method.  Test  sequence  is  ARMA4  CL.  (a) 
Model  order  (4,4)  with  15  dB  of  added  noise  and  (b)  model  order  (8,8) 
with  5  dB  of  added  noise. 
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Figure  3.19:  Model  pole  scatter  plots  of  twenty  trials  on  test  sequence 
ARMA4  CL  illustrating  (a)  the  ability  of  correlation  domain  prefiltering 
to  improve  the  resolution  of  an  LSMYWE  initial  estimate,  order  (8,8),  5 
dB  added  noise  and  (b)  the  ability  of  Akaike  MLE  to  improve  the  resolution 
of  a  Prony's  method  initial  estimate,  model  order  (6,2),  30  dB  added  noise. 
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Figure  3.20:  Pole-zero  model  of  the  ARMA4  CL  test  sequence  with  5  dB 
of  additive  white  noise  using  LSMYWE  and  the  smallest  singular  value 
discarded.  The  MA  part  was  found  using  least  squares  identification. 
This  example  corresponds  to  the  example  in  Figure  (3.15a).  (a)  The  noisy 
sequence  and  (b)  the  estimated  model  impulse  response  and  the  original 
noiseless  sequence.  Model  order  (8,8). 
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IV.  LINEAR  MODELING  OF  ACOUSTIC 

TRANSIENTS 

A.  ACOUSTIC  DATA— GENERAL 

The  laboratory  generated  acoustic  transient  data  available  for  this  study  were 
generated  in  six  different  ways  using  ordinary  laboratory  objects.  The  data  records 
and  their  method  of  generation  are  summarized  in  Table  4.1.  Each  transient  will 
be  refered  to  by  the  object  used  or  the  action  performed  to  generate  that  particular 
transient  such  as  'book'  or  'slam'.  The  six  data  records  modeled  in  this  section  are 
shown  in  Figures  4.1a-f.  The  range  of  data  that  was  actually  modeled  is  listed  in 
Table  4.1.  The  sampling  rate  for  the  data  is  unknown  but  is  not  required  since  there 
is  no  need  to  infer  the  specific  characterics  of  the  acoustic  source  or  the  acoustic 
enviroment. 


Data  Name 

Generation  Technique 

Indices  of  Modeled  Data 
(beginrend) 

Book 

Dropped  Book 

(55:454) 

Slam 

Slammed  Book 

(15:414) 

Plate 

Struck  Metal  Plate 

(5:704) 

Ruler 

Book  Struck  with  Ruler 

(501:650) 

Wood 

Clapped  Wooden  Blocks 

(171:320) 

Wrench 

Dropped  Wrench 

(101:250) 

TABLE  4.1:  Summary  of  Acoustic  Transient  Data  and  method  of  genera- 
tion. 


B.     ACOUSTIC  TRANSIENT  MODELING  RESULTS 

Iterative  prefiltering  and  correlation  domain  iterative  prefiltering  were  found  to 
yield  the  most  effective  time  domain  match  (in  terms  of  squared  error)  between  the 
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original  signal  and  the  impulse  response  of  the  estimated  model.  LSMYWE  with 
the  smallest  singular  value  set  to  zero,  as  in  Figures  3.13a,b,  was  used  to  find  the 
initializing  model  poles  for  the  iterative  algorithms.  Removing  the  smallest  singular 
value  allowed  the  LSMYWE  algorithm  to  place  zeros  much  closer  to  the  unit  circle 
than  was  possible  when  all  singular  values  were  used.  Prony's  method  and  the  Akaike 
MLE  algorithm  were  also  used  to  model  laboratory  transients.  However,  their  perfor- 
mance was  generally  poor  and  they  will  not  appear  in  the  remainder  of  this  chapter. 
As  noted  in  Chapter  III,  the  problem  with  Prony's  method  is  that  it  cannot  easily 
model  highly  resonant  frequencies,  that  is,  zeros  close  to  the  unit  circle,  in  the  pres- 
ence of  even  small  amounts  of  noise  unless  a  large  number  of  excess  poles  are  used 
and  numerous  singular  values  are  discarded.  This  procedure  is  burdensome  when 
initializing  iterative  algorithms  which  do  not  require  these  excess  poles.  The  primary 
difficulty  with  Akaike  MLE  is  its  inability  to  handle  zeros  outside  the  unit  circle.  The 
initializing  model  zeros  were  found  using  Shank's  method  which  was  by  far  the  most 
robust  algorithm  for  finding  numerator  coefficients  of  the  methods  tested  in  Chapter 
Three.  In  addition  to  the  time  domain  plots  of  model  impulse  response,  a  pole-zero 
model  spectrum  and  pole-zero  plot  is  be  provided  for  each  transient  so  that  their 
differing  characteristics  can  be  observed.  All  spectra  were  generated  by  squaring  the 
FFT  magnitude  of  either  the  model  coefficients  or  the  signal  being  modeled.  Model 
order  was  chosen  based  on  the  best  educated  guess  of  the  author  in  accordance  with 
the  recommendations  in  Chapter  III,  and  augmented  by  trial  and  error. 

The  LSMYWE  model  used  to  initialize  the  iterative  preflltering  algorithm  for 
the  Slam  transient  is  shown  in  Figures  4.2a  and  4.3a.  The  corresponding  iterative 
preflltering  model  shown  in  Figures  4.2b  and  4.3b.  Figure  4.4a  shows  the  iterative 
preflltering  model  spectrum.  Figure  4.4b  illustrates  a  characteristic  of  the  itera- 
tive preflltering  algorithm  that  was  observed  throughout  the  modeling  of  acoustic 
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transients.  Namely,  when  an  excess  number  of  model  parameters  are  used,  the  error 
between  the  model  impulse  response  and  the  signal  being  modeled  increases  dramat- 
ically, mainly  at  the  beginning  of  the  signal.  The  subsequent  application  of  Shank's 
method  is  effective  in  reducing  this  initial  error.  This  effect  is  shown  in  Figure  4.5a. 
In  general,  the  application  of  Shank's  method  as  the  last  step  in  the  modeling  process 
was  found  to  reduce  the  the  mean  squared  error  between  the  original  signal  and  the 
model  impulse  response  to  some  degree  for  all  modeling  methods.  The  sensitivity  to 
excess  parameters  shown  by  iterative  prefiltering  does  not  affect  correlation  domain 
iterative  prefiltering.  Note  that  excess  model  zeros  cause  no  difficulties  in  Figure 
4.5b.  For  the  Book  transient,  two  modeling  trials  are  shown.  Figures  4.6  and  4.7 
show  the  best  time  domain  match  obtained  using  iterative  prefiltering  and  also  an 
LSMYWE,  Shank's  method,  respectively.  Although  the  LSMYWE  model  does  not 
achieve  as  effective  an  impulse  response  match  as  iterative  prefiltering,  with  SVD  it 
is  more  sensitive  to  the  low  energy,  high  frequency  component  present  in  the  Book 
transient  at  approximately  |.  The  best  model  of  Ruler  is  shown  in  Figure  4.8.  The 
model  spectra  for  Book  and  Ruler  appear  in  Figure  4.9. 

The  assumption  that  the  signal  being  modeled  is  a  system  impulse  response 
is  problematic  for  the  Plate,  Wood,  and  Wrench  signals  since  they  do  not  exhibit 
the  rapid  decay  usually  associated  with  an  impulse  response.  However,  it  is  still 
possible  to  achieve  a  resonable  time  domain  match  over  small  segments  of  each  signal. 
This  result  is  illustrated  in  the  models  of  in  Figures  4.10,  4.11,  4.13,  and  4.14.  The 
model  spectra  for  Wood  and  Wrench  are  shown  in  Figure  4.15.  The  reason  that 
correlation  domain  iterative  prefiltering  was  used  for  the  Plate,  Wood,  and  Wrench 
signals  instead  of  standard  iterative  prefiltering  can  be  illustrated  by  comparing  the 
two  techniques  on  a  segment  of  the  Plate  signal.  Although  the  impulse  response  error 
of  the  two  models  in  Figures  4.10  and  4.11  are  nearly  identical,  Figure  4.12  shows 
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that  correlation  domain  iterative  prefiltering  clearly  outperformed  standard  iterative 
prefiltering  in  reproducing  the  spectrum  of  the  original  sequence.  In  fact,  at  the  more 
suitable  model  order  of  (16,16),  iterative  prefiltering  would  not  converge  but  instead 
oscillated  in  a  region  of  convergence  for  any  model  order  over  (12,12).  Figure  4.16b 
shows  the  model  impulse  response  obtained  when  the  large  segment  of  Plate  shown 
in  Figure  4.16a  is  modeled  using  iterative  prefiltering.  Although  the  time  domain  and 
spectral  properties  of  the  model  relative  to  the  original  signal  are  considerable  poorer 
than  those  obtained  for  a  short  segment,  the  model  does  clearly  share  many  of  the 
features  of  the  original  signal. 

C.     ACOUSTIC  SIGNAL  MODELING  SUMMARY 

The  modeling  results  obtained  in  the  previous  section  of  this  Chapter  indicate 
the  possible  utility  of  pole-zero  modeling  algorithms  with  regard  to  modeling  tran- 
sient signals.  Signals  with  decaying  narrowband  components  (e.g.  Slam,  Book,  and 
Ruler)  and  signals  with  substained  narrowband  components  (e.g.  Plate,  Wood,  and 
Wrench)  can  be  modeled  as  the  impulse  response  of  a  rational  linear  system.  Robust 
modeling  algorithms  are  available  which  can  effectively  deal  with  the  many  uncertain- 
ties associated  with  real  world  signals.  Although  the  goal  of  achieving  an  exact  time 
domain  match  between  the  original  signal  and  the  pole-zero  model  impulse  response 
was  not  realized  for  any  of  the  acoustic  signals  in  this  chapter,  in  all  cases  the  degree 
of  match  obtained  clearly  indicates  that  many  signal  characteristics  are  described  by 
the  model. 
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Figure  4.1:   Laboratory  generated  acoustic  transient  data:    (a)  Slam  and 
(b)  Book. 
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Figure  4.1:   continued  Laboratory  generated  acoustic  transient  data:    (c) 
Ruler  and  (d)  Plate. 
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Figure  4.1:   continued  Laboratory  generated  acoustic  transient  data:    (e) 
Wood  and  (f)  Wrench. 
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Figure  4.2:  Modeling  the  transient  Slam — model  impulse  response  versus 
the  original  signal.  Model  order  (6,8).  (a)  LSMYWE  and  Shank's  method 
and  (b)  iterative  prefiltering  initialized  with  (a). 
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Figure  4.3:  Modeling  the  transient  Slam — model  pole-zero  plots.  Model 
order  (6,8).  (a)  LSMYWE  with  SVD  and  Shank's  method  and  (b)  iterative 
prefiltering  initialized  with  (a). 


73 


xlO' 


2000 


1000  - 


12        3 
frequency  (rad) 

original  signal 

model  impulse  response 


-1000 


-2000 


100 


200 


n 


300 


400 


(a) 


(b) 


Figure  4.4:  Modeling  the  transient  Slam — model  spectrum  and  and  over- 
parametrized  iterative  prefiltering  model  impulse  response  versus  the  orig- 
inal signal,  (a)  Model  spectrum  for  iterative  prefiltering  with  model  order 
(6,8)  and  (b)  iterative  prefiltering  for  model  order(6,12).  Note  how  excess 
model  parameters  cause  error  at  the  beginning  of  the  iterative  prefiltering 
model. 
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Figure  4.5:  Modeling  the  transient  Slam — model  impulse  response  versus 
the  original  signal,  overparameterized  modeling  effects,  (a)  The  applica- 
tion of  Shanks  method  to  the  iterative  prefiltering  model  of  order  (6,12) 
reduces  the  error  at  the  beginning  of  the  model  impulse  response  and 
(b)  correlation  domain  iterative  prefiltering  for  model  order  (6,12)  can 
effectively  use  the  excess  model  parameters. 
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Figure  4.6:  Modeling  the  transient  Book — model  impulse  response  versus 
the  original  signal.  Model  order  (6,6).  (a)  LSMYWE  with  SVD  and 
Shank's  method  and  (b)  iterative  prefiltering  initialized  with  (a).  Note 
the  sensitivity  of  LSMYWE  with  SVD  to  the  high  frequency  components 
present. 
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Figure  4.7:  Modeling  the  transient  Book — model  pole-zero  plots.  Model 
order  (6,6).  (a)  LSMYWE  with  SVD  and  Shank's  method  and  (b)  iterative 
prefiltering  initialized  with  (a).  Note  the  sensitivity  of  LSMYWE  with 
SVD  to  the  high  frequency  pole  pair. 
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Figure  4.8:  Modeling  the  transient  Ruler — Model  order  (6,16).  (a)  Iter- 
ative prefiltering  model  impulse  response  versus  the  original  signal  and 
(b)  the  corresponding  model  pole-zero  plot.  Note  that  the  two  beating 
sinusoids  have  been  effectively  modeled. 
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Figure  4.9:  Modeling  the  transient  Book  and  Ruler — model  spectrums. 
(a)  The  iterative  prefiltering  model  of  order  (6,6)  for  the  transient  Book 
and  (b)  the  iterative  prefiltering  model  of  order(6,12)  for  the  transient 
Ruler. 
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Figure  4.10:  Modeling  the  transient  Plate — Model  order  (12,12).  (a)  It- 
erative prefiltering  model  impulse  response  versus  the  original  signal  and 
(b)  the  corresponding  model  pole-zero  plot. 
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Figure  4.11:  Modeling  the  transient  Plate — Model  order  (16,16).  (a)  Cor- 
relation domain  iterative  prefiltering  model  impulse  response  versus  the 
original  signal  and  (b)  the  corresponding  model  pole-zero  plot. 
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Figure  4.12:  Modeling  the  transient  Plate — model  spectrums  versus  the 
original  signal  spectrums.  (a)  The  iterative  prefiltering  model  of  order 
(12,12)  and  (b)  the  correlation  domain  iterative  prefiltering  model  of  order 
(16,16) 
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Figure  4.13:  Modeling  the  transient  Wood — Model  order  (16,16).  (a) 
Correlation  domain  iterative  prefiltering  model  impulse  response  versus 
the  original  signal  and  (b)  the  corresponding  model  pole-zero  plot. 
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Figure  4.14:  Modeling  the  transient  Wrench — Model  order  (8,8).  (a)  Cor- 
relation domain  iterative  prefiltering  model  impulse  response  versus  the 
original  signal  and  (b)  the  corresponding  model  pole-zero  plot. 
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Figure  4.15:  Modeling  the  transient  Wood  and  Wrench — model  spectrums. 
(a)  The  iterative  prefiltering  model  of  order  (16,16)  for  the  transient  Wood 
and  (b)  the  iterative  prefiltering  model  of  order(8,8)  for  the  transient 
Wrench. 
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inal Plate  segment  and  (b)  the  model  impulse  response. 
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V.  CONCLUSIONS 

A.     PERFORMANCE  COMPARISON  SUMMARY 

The  modeling  of  signals  as  the  impulse  response  of  a  linear  pole- zero  system  is 
an  important  tool  in  digital  signal  processing.  One  natural  area  for  applying  such 
an  approach  is  in  the  modeling  of  transient,  impulse  response-like  waveforms.  The 
specific  approach  taken  in  this  thesis  was  to  determine  which  pole-zero  modeling  al- 
gorithms are  most  suited  to  modeling  complex,  real  world  transient  waveforms.  The 
modeling  criterion  emphasized  is  to  obtain  the  best  (least  squares  error)  time  domain 
match  between  the  model  impulse  response  and  the  original  signal.  This  criterion 
was  adopted  because  it  provides  a  degree  of  signal  characterization  that  is  a  step  be- 
yond normal  power  spectrum  estimation.  Indeed,  the  strength  of  pole- zero  models  is 
their  ability  to  describe  not  only  the  resonances  present  (model  poles),  but  also  how 
these  resonances  are  related  (model  zeros).  Because  of  the  widely  varying  character- 
istics anticipated  for  real  world  signals,  a  key  evaluation  criterion  for  any  modeling 
algorithm  is  robustness.  In  particular,  algorithms  must  be  effective  in  the  presence  of 
signal  degrading  effects  like  noise  and  model  degrading  effects  such  as  unknown  model 
order.  The  performance  of  several  selected  algorithms  were  compared  for  known  im- 
pulse response  test  sequences  in  Chapter  III.  The  modeling  experience  gained  in  these 
experiments  was  then  applied  to  modeling  laboratory  generated  acoustic  transient 
datain  in  Chapter  IV  as  a  test  of  'real  world'  effectiveness. 

Four  basic  algorithms  were  chosen  for  comparison:  Prony's  method,  the  least 
squares  modified  Yule- Walker  equations  (LSMYWE),  iterative  prefiltering,  and  the 
Akaike  maximum  likelihood  estimator  (MLE).  An  algorithm  which  is  an  extension  of 
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iterative  prefiltering  into  the  correlation  domain  was  also  presented.  For  those  meth- 
ods in  which  the  model  poles  and  zeros  are  determined  separately  (Prony's  method 
and  LSMYWE),  four  methods  for  determining  model  zeros  (i.e.,  transfer  function  nu- 
merator coefficients)  were  considered:  the  upper  partition  of  Prony's  method,  spectral 
factorization,  Durbin's  method,  and  Shank's  method.  The  major  conclusions  of  the 
algorithm  performance  comparisons  conducted  in  Chapter  III  and  Chapter  IV  are  as 
follows: 

1.  Algorithms  that  are  unable  to  model  zeros  outside  the  unit  circle  (Durbin's 
method,  Akaike  MLE)  have  limited  versatility  when  modeling  arbitrary  tran- 
sient waveforms.  All  the  acoustic  transients  in  Chapter  IV  required  a  non- 
minimum  phase  model  to  obtain  the  best  time  domain  match. 

2.  The  most  robust  and  effective  method  for  finding  zeros  that  gives  the  best  least 
squares  time  domain  match  was  found  to  be  Shank's  method.  In  fact,  applying 
Shank's  method  as  a  last  step  improved  the  final  model  of  all  algorithms  to 
some  degree.  The  upper  partition  of  Prony's  method  and  spectral  factorization 
are  not  very  useful  because  of  their  extreme  sensitivity  to  noisy  or  time  shifted 
signals. 

3.  Prony's  method  and  Akaike  MLE  have  difficulty  modeling  signals  in  which 
additive  noise  is  present.  Even  small  amounts  of  additive  noise  causes  a  dramatic 
loss  of  pole  frequency  resolution  for  Prony's  method.  The  use  of  excess  poles 
and  singular  value  decomposition  were  found  to  be  effective  in  overcoming  these 
effects  but  these  methods  depend  on  a  knowledge  of  correct  model  order.  Also, 
unlike  spectrum  estimation,  excess  poles  must  be  retained  for  time  domain 
matching.  This  is  in  direct  contrast  to  the  parsimonious  use  of  model  parameters 
normally  provided  by  a  pole-zero  model.  The  difficulties  encountered  with  the 
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Akaike  MLE  were  computational  in  nature;  for  noisy  signals  the  Akaike  MLE 
algorithm  usually  terminated  prematurely  when  either  a  non-minimum  phase 
model  was  encountered  during  an  iteration  or  when  a  numerical  overflow  was 
induced  by  an  unstable  model  estimate. 

4.  LSMYWE  with  singular  value  decomposition  and  iterative  prefiltering  (includ- 
ing the  correlation  domain  version)  were  found  to  be  the  most  effective  algo- 
rithms for  modeling  a  signal  when  additive  noise  is  present.  If  the  true  model 
order  of  a  system  is  unknown,  it  is  best  to  discard  only  one  singular  value. 
Almost  all  the  resolution  gain  occurs  with  the  first  singular  value.  Discarding 
additional  singular  values  is  intended  to  reduce  the  variance  of  excess  poles  but 
it  will  cause  poles  to  be  biased  if  all  model  poles  are  necessary  for  signal  model- 
ing. Both  of  these  methods  demonstrated  the  consistent  ability  to  model  poles 
very  close  to  the  unit  circle.  This  capability  was  essential  when  modeling  the 
acoustic  transients  used  in  this  thesis. 

Combining  the  above  observations  leads  to  our  recommended  strategy  for  mod- 
eling an  arbitrary  transient  signal.  First,  use  LSMYWE  with  one  singular  value  re- 
moved and  Shank's  method  to  find  an  initial  model  estimate.  Next,  use  this  estimate 
to  initialize  either  iterative  prefiltering  or  correlation  domain  iterative  prefiltering.  Fi- 
nally, if  desired,  apply  Shank's  method  to  optimize  the  time  domain  fit  of  the  model 
zeros. 

B.     RECOMMENDATIONS  FOR  FUTURE  STUDY 

The  results  of  Chapter  IV  indicate  that  there  are  pole-zero  modeling  algorithms 
available  that  are  sufficiently  robust  to  be  useful  for  modeling  many  complex,  real 
world  transients.  A  number  of  issues  regarding  the  pole-zero  modeling  of  transient 
signals  and  the  application  of  such  models  require  further  study.  These  issues  include: 
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1.  The  degree  to  which  the  noise  performance  of  correlation  domain  iterative  pre- 
filtering  may  be  improved  by  introducing  noise  compensation  of  the  autocorre- 
lation sequence  requires  study. 

2.  A  study  of  the  effectiveness  of  model  order  determination  techniques  when  ap- 
plied to  transient  modeling  would  facilitate  more  effective  use  of  transient  mod- 
eling techniques. 

3.  The  effectiveness  of  iterative  prefiltering  and  correlation  domain  prefiltering  as 
a  spectral  estimation  technique  should  be  explored  further. 

4.  The  ability  of  pole-zero  models  to  describe  the  time  domain  characteristics  of 
a  signal  could  aid  in  the  detection  of  signals.  Such  an  application  needs  to  be 
persued. 

5.  A  study  of  the  structural  relationship  of  pole- zero  models  to  the  specific  systems 
that  generate  transients  may  increase  the  usefulness  of  such  models. 
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